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4.  INTRODUCTION 

The  HBCU/MI  Partnership  Training  Award  project,  “ Photonic  Breast  Tomography  and 
Tumor  Aggressiveness  Assessment is  designed  to  establish  a  breast  cancer  training  and  research 
program  at  the  City  College  of  New  York  (CCNY)  through  close  collaboration  with  the 
researchers  at  the  Memorial  Sloan  Kettering  Cancer  Center  (MSKCC).  The  focus  of  the  training 
component  of  the  project  is  to  familiarize  the  CCNY  researchers  who  happen  to  be  physical 
scientists  and  engineers  to  the  biological  aspects  of  cancer  research  through  attending  relevant 
courses,  and  cancer  research  practicum  through  laboratory  rotations.  The  objectives  of  the 
research  component  of  the  project  are  to  develop  optical  imaging  and  spectroscopic  approaches 
to  (a)  distinguish  between  aggressive  and  slow  growing,  metastatic  and  non-metastatic  tumors, 
(b)  non-invasively  detect  and  diagnose  breast  tumors  at  early  stages  of  growth. 

During  the  fourth  reporting  period  (June  15,  2010  -  June  14,  2011)  covered  by  this  report, 
the  major  research  thrust  was  on  developing  near- infrared  (NIR)  light-based  experimental 
methods  and  numerical  algorithms  for  detection  of  targets  in  breast  tissue  simulating  model 
medium  with  the  goal  of  detecting  breast  tumors  at  early  stages  of  growth.  The  trainees  also 
continued  learning  the  magnetic  resonance  spectroscopic  imaging  (MRSI)  with  selective 
multiple-quantum  coherence  transfer  (SelMQC)  sequence  to  explore  its  potential  to  investigate 
the  glycolytic  activity  of  tumors. 

5.  BODY 

We  have  made  substantial  progress  in  developing  non-invasive  near-infrared  optical 
imaging  modalities  for  early  detection  of  breast  cancer  ( Specific  Aim  4)  that  we  started  during 
Year  1  of  the  project  [1]  and  pursued  and  reported  in  the  subsequent  years  [2,3],  The  goal  of  the 
research  is  to  develop  optical  spectroscopy  and  imaging  approaches  that  use  the  near-infrared 
light  to  obtain  three-dimensional  (3-D)  tomographic  images  of  human  breast  to  enable  detection, 
localization,  and  possible  diagnosis  of  tumor(s)  in  the  breast.  Early  detection  of  breast  tumor 
with  requisite  sensitivity  and  specificity  is  a  daunting  task  and  we  continued  development  of 
different  approaches,  which  include: 

•  Diffuse  Optical  Imaging  using  Decomposition  Methods; 

•  Time  Reversal  Optical  Tomography  (TROT);  and 

•  Finite  element  method  (FEM). 

We  provide  a  brief  outline  of  activities  and  accomplishments  in  these  areas,  and  refer  to 
appended  materials  for  detailed  description  where  applicable. 

5.1.  Diffuse  Optical  Imaging  Using  Decomposition  Methods 

We  have  previously  reported  on  the  development  of  Optical  Tomography  using  Independent 
Component  Analysis  (OPTICA),  extended  it  to  include  multi-wavelength  probing  and  explored 
its  efficacy  on  a  more  realistic  model  breast  ( Specific  Aim  4,  Task  #15  and  Task  #16)  [2-4], 
Diffuse  optical  imaging  (DOI)  for  detection  and  retrieval  of  location  information  of  targets  in  a 
highly  scattering  turbid  medium  may  be  treated  as  a  “blind  source  separation”  (BSS)  problem 
[5],  Independent  Component  Analysis  (ICA),  the  basis  for  OPTICA,  is  one  of  the  different 
matrix  decomposition  methods  for  solving  the  BSS  problem  and  retrieving  desired  information. 
Other  decomposition  methods  include  Principal  Component  Analysis  (PCA)  [6]  and  Non- 
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negative  Matrix  Factorization  (NMF)  [7].  The  three  algorithms  have  different  assumptions, 
which  may  lead  to  different  favored  conditions. 

ICA  assumes  the  signals  from  different  targets  to  be  independent  of  each  other,  and 
optimizes  a  relevant  measure  of  independence  to  obtain  the  ICs  associated  with  different  targets 
[4],  The  position  co-ordinates  of  targets  in  three  dimensions  are  determined  from  the  individual 
components  separately.  PCA  [6]  assumes  that  the  PCs  contributing  to  the  signal  are  uncorrelated 
and  explain  the  most  variance  in  the  signal.  NMF  [7]  seeks  to  factorize  a  matrix  into  two  non¬ 
negative  matrices  (component  signals  and  weights)  and  requires  the  contributions  to  signal  and 
the  weights  of  the  components  to  be  non-negative.  It  does  not  imply  any  relationship  between  the 
components.  Our  objective  was  to  test  and  compare  the  efficacy  of  these  three  approaches  in 
solving  the  DOI  problem.  We  used  both  simulative  data  and  experimental  data  for  absorptive  and 
scattering  targets  embedded  in  model  scattering  media.  Details  of  the  theoretical  formalism, 
numerical  algorithms,  simulation  and  experiment  are  detailed  in  the  pre-print  of  a  submitted 
paper  {attached  as  Appendix  3);  only  a  brief  overview  is  presented  below. 

Blind  source  separation  (BSS),  also  known  as  blind  signal  separation  is  a  general  problem  in 
information  theory  that  seeks  to  separate  the  contributions  from  different  sources  to  the  measured 
signal,  which  is  a  weighted  mixture  of  signals  from  those  sources.  Assuming  the  source  signals 
are  linearly  mixed,  the  BSS  problem  can  be  presented  in  matrix  notation  as,  X=AS,  where  X 
represent  measured  signal,  A  is  a  mixing  or  weighting  matrix,  and  S  represents  signals  from  the 
sources.  The  objective  of  BSS  is  to  retrieve  the  source  signals  and  their  weights  from  the 
measured  signal.  Due  to  the  lack  of  prior  knowledge  of  the  source  signals,  statistical  analysis 
methods,  such  as  ICA,  PCA  and  NMF  are  used  to  retrieve  source  information. 

In  DOI  one  measures  the  signal  at  the  sample  boundary,  which  is  a  weighted  mixture  of 
contributions  from  embedded  targets.  The  signal  is  the  perturbation  in  the  light  intensity 
distribution  at  the  sample  boundary,  and  ideally  is  the  difference  between  the  signal  recorded 
with  the  target(s)  and  that  without  the  targets.  In  practice  signal  without  the  targets  is  reasonably 
approximated  by  an  average  of  the  signals  recorded  for  different  scanning  positions  of  the 
sample.  DOI  is  thus  a  BSS  problem  in  the  optical  domain,  and  the  problem  was  investigated 
using  three  matrix  decomposition  approaches,  ICA,  PCA  and  NMF,  for  absorptive  and  scattering 
targets. 

In  simulation,  the  sample  was  considered  to  be  a  40-mm  thick  uniform  slab  of  scattering 
medium,  as  shown  in  Fig.  1  of  Appendix  3.  Its  absorption  and  diffusion  coefficients  were  taken 
to  be  pa  =  0.003  mm'1  and  D  =  1/3  mm,  respectively,  which  are  similar  to  the  average  value  of 
those  parameters  for  human  breast  tissue.  An  absorptive  and  a  scattering  point  target  were  placed 
at  (50,  60,  15)  mm  and  (30,  30,  25)  mm,  respectively.  The  absorption  coefficient  of  the 
absorptive  target  was  set  to  be  higher  than  the  background  with  A pa  =  0.01  mm'1,  while  the 
diffusion  coefficient  was  taken  to  be  the  same  as  that  of  background.  The  diffusion  coefficient  of 
the  scattering  target  was  set  to  be  lower  than  the  background  (higher  scattering  coefficient)  with 
AD  =  -1  mm,  while  the  absorption  coefficient  was  taken  to  be  the  same  as  the  background.  The 
incident  CW  beam  was  assumed  to  step  scan  the  sample  at  21x21  grid  points  covering  80x80 
mm2  area,  with  a  step  size  of  4  mm.  Light  on  the  opposite  side  was  recorded  at  41  x41  grid  points 
covering  the  same  area.  Multiplicative  Gaussian  noise  of  5%  was  added  to  the  simulated  data. 
The  data  matrix  X  was  then  obtained,  and  analyzed  using  the  three  different  algorithms. 
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To  evaluate  the  efficacy  of  the  decomposition  methods  in  experiments,  we  carried  out  two 
different  experiments  with  two  different  samples.  The  first  sample  used  a  250  mm  x  250  mm  x 
50  mm  transparent  plastic  container  filled  with  Intralipid- 10%  suspension  in  water  as  the 
background  medium.  The  concentration  of  Intralipid- 10%  was  adjusted  to  provide  [8,9]  an 
absorption  coefficient  of  pa  ~  0.003  mm-1,  and  a  transport  mean  free  path  lt  ~  1.43  mm  at  785 
nm.  The  second  sample  used  a  similar  container  with  dimension  of  250  mm  x  250  mm  x  60  mm 
filled  with  Intralipid-20%  suspension  in  water.  The  concentration  of  Intralipid-20%  was  adjusted 
to  provide  [8,  9]  pa  ~  0.003  mm'1,  and  lt  ~  1  mm  at  790  nm.  These  optical  parameters  of  the 
medium  were  selected  to  be  similar  to  those  for  human  breast  tissue.  The  thickness  of  the 
samples  was  also  comparable  to  that  of  a  typical  compressed  female  human  breast. 

In  the  first  experiment,  two  absorptive  targets  were  embedded  in  the  medium.  The  targets 
were  ~  10-mm  diameter  glass  spheres  filled  with  a  solution  of  Indo-cyanine  green  dye  in  water. 
The  absorption  coefficient  pa  was  adjusted  to  be  1.15  mm'1  at  785  nm,  with  jus  approximately  the 
same  as  that  of  the  background  medium.  The  targets  were  placed  at  (57.2,  18.1,  20.0)  mm  and 
(19.9,  48.1,  25.0)  mm,  respectively.  In  the  second  experiment,  two  scattering  targets  were 
embedded,  which  were  also  ~  10  mm  diameter  glass  spheres,  filled  with  Intralipid-20% 
suspension  in  water.  The  transport  mean  free  path,  lt  was  adjusted  to  be  0.25  mm,  with  scattering 
coefficient  ps  ~  11  mm'1,  and  absorption  coefficient  pa  same  as  the  background  medium.  The 
targets  were  placed  in  the  mid-plane  (z  =  30  mm)  in  the  container  with  a  lateral  distance  of  40 
mm  from  each  other  (center  to  center). 

The  experimental  setup  (shown  schematically  in  Fig.  6  of  Appendix  3),  was  based  on  what 
we  assembled  {Specific  Aim  4,  Task  #14)  earlier  in  the  project  [1].  A  10-mW  785-nm  diode  laser 
beam  was  used  to  illuminate  the  first  sample,  while  a  100-mW  790-nm  diode  laser  beam  was 
used  for  the  second  sample.  The  input  surface  (source  plane)  of  the  samples  was  scanned  across 
the  laser  beam  in  an  x-y  array  of  grid  points  to  realize  the  multi-source  interrogation  of  the 
samples.  The  transmitted  light  from  the  exit  surface  (detector  plane)  was  recorded  by  a  1024 
pixel  xl024  pixel  (pixel  size  =  24  pm)  CCD  camera  (Photometries  CH350)  equipped  with  a  60- 
mm  focal-length  camera  lens.  Each  pixel  of  the  CCD  camera  can  be  considered  to  be  a  detector 
implementing  the  multi-detector  signal  acquisition  arrangement.  A  set  of  16-bit  1024  pixel 
xl024  pixel  images  were  acquired.  The  two  samples  were  scanned  in  an  array  of  11x12  and 
11x15  grid  points,  respectively,  with  a  step  size  of  5  mm  in  both  cases.  The  processes  of 
scanning  and  data  acquisition  were  controlled  by  a  personal  computer.  At  all  scan  positions,  raw 
transillumination  images  of  the  samples  were  recorded  by  the  computer  for  further  analysis.  A 
16-bit  1024  x  1024  pixels  CCD  camera  recorded  the  transmitted  signal.  The  sample  was  scanned 
across  the  laser  beam  in  a  16x26  x-y  array  of  grid  points  and  a  two-dimensional  transmission 
image  was  recorded  for  each  position  for  each  wavelength  to  meet  the  multi-source  multi¬ 
detector  (each  pixel  of  a  CCD  camera  being  a  detector  element)  imaging  arrangement.  The 
resulting  data  was  analyzed  using  the  numerical  algorithm  of  ICA,  PCA  and  NMF.  The  details  of 
the  experimental  arrangement,  data  acquisition,  processing  and  analysis  methods  are  presented  in 
attached  Appendix  3. 

The  key  result  from  simulated  data  is  that  the  positions  and  optical  strengths  of  the  targets 
retrieved  by  ICA,  PCA  and  NMF  algorithms  in  excellent  agreement  with  the  known  values. 

The  results  from  the  experiments  using  absorptive  and  scattering  targets  demonstrate  that 
overall,  all  three  algorithms  detect  and  locate  the  scattering  and  the  absorptive  targets  with  good 
accuracy,  the  maximum  deviation  of  any  one  coordinate  from  the  known  value  being  -3  mm. 
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The  (x,  y,  z)  positions  of  the  targets  were  retrieved  with  good  accuracy.  The  decomposition 
provided  by  ICA  is  “cleaner”  than  that  of  the  PCA.  PCA  did  not  clearly  separate  the  two 
absorptive  targets  used  in  the  first  experiment.  NMF  decomposition  seems  to  provide  residue- 
free  “cleaner”  images  than  the  other  two  methods  in  this  study.  Appendix  3  provides  further 
discussion  on  implication  of  results. 

5.2.  Time  Reversal  Optical  Tomography 

We  have  been  pursuing  development  of  the  Time  Reversal  Optical  Tomography  (TROT) 
approach,  [2,  3]  in  our  on-going  quest  for  fast  and  accurate  methods  for  detection  and  localization 
of  tumours  in  breast,  and  for  detection  of  margins  during  surgical  removal  of  breast  tumours 
(Specific  Aim  4). 

The  time  reversal  (TR)  invariance,  the  basic  symmetry  that  commonly  holds  in  microscopic 
physics,  forms  the  basis  for  macroscopic  imaging  [10]  in  TROT.  TROT  also  adapts  and 
incorporates,  in  optical  domain,  the  vector  subspace  classification  method,  Multiple  Signal 
Classification  (MUSIC).  MUSIC  was  developed  by  Devaney  and  co-workers  for  finding  the 
location  of  scattering  targets  whose  size  is  smaller  than  the  wavelength  of  acoustic  waves  or 
electromagnetic  waves  (radar)  used  for  probing  the  homogeneous  or  inhomogeneous  background 
medium  in  which  the  targets  were  embedded  [11,  12].  In  optical  imaging  application,  a  response 
matrix  represents  the  transport  of  light  from  multiple  sources  through  a  turbid  medium  with 
embedded  targets  to  an  array  of  detectors.  The  response  matrix  is  constructed  from  the 
experimental  data  that  represent  the  perturbation  in  light  intensity  distribution  on  the  sample 
boundary  due  to  the  presence  of  target(s).  The  ‘time  reversal  (TR)  matrix’  is  constructed  by 
multiplying  the  response  matrix  by  its  transpose  matrix  for  continuous  wave  illumination  (by 
adjoint  matrix  for  frequency-domain  case).  Mathematically,  the  TR  matrix  is  equivalent  to  transfer 
of  light  from  multiple  sources  through  a  turbid  medium  with  embedded  targets  to  an  array  of 
detectors,  and  back,  and  is  similar  to  the  time-reversal  matrix  used  in  the  general  area  of  array 
processing  for  acoustic  and  radar  time-reversal  imaging  [12].  The  eigenvalue  equation  of  TR 
matrix  is  solved,  and  the  signal  and  noise  are  separated  into  orthogonal  subspaces,  using  an  L- 
curve  regularization  method.  Then  a  pseudo  spectrum  is  calculated  directly  for  all  voxels  in  the 
sample  using  MUSIC  [11,  12].  Pseudo  tomographic  images  can  be  generated  using  pseudo  values. 
Locations  of  targets  are  determined  by  the  global  maxima  (or  local  maxima  in  low  SNR  case) 
components  in  the  pseudo  spectrum.  Details  of  the  theoretical  formalism  and  numerical  algorithms 
of  TROT  along  with  tests  of  its  efficacy  using  simulative  and  experimental  data  have  been  detailed 
in  a  manuscript  submitted  for  publication  in  and  is  attached  to  this  report  as  Appendix  4.  We  will 
present  only  the  key  results  in  the  body  of  this  report,  and  refer  to  Appendix  4  for  further  details. 

First  we  tested  the  potential  of  TROT  in  an  ideal  situation  using  a  challenging  problem  in 
simulation  that  involved  detecting  and  locating  6  targets  embedded  in  a  40-mm  thick  breast  tissue- 
simulating  scattering  medium.  The  six  point-like  absorptive  targets,  with  absorption  coefficient 
difference  of  A pa  =  0.01  mm-1  from  the  background,  were  placed  at  A  (24  mm,  26  mm,  9  mm),  B 
(38  mm,  38  mm,  15  mm),  C  (38  mm,  38  mm,  21  mm),  D  (40  mm,  38  mm,  21  mm),  E  (44  mm,  42 
mm,  21  mm)  and  F  (30mm,  30mm,  31  mm),  respectively.  The  origin  (0  mm,  0  mm,  0  mm)  was 
located  at  the  upper-left  comer  of  the  input  boundary  (source  boundary)  of  the  sample.  As  can  be 
seen  from  the  assigned  coordinates,  targets  C  and  D  are  located  at  two  adjacent  grid  points,  and  are 
close  to  target  E,  and  these  three  targets  are  located  in  the  same  z  layer.  Consequently,  targets  C 
and  D  are  expected  to  be  very  difficult  to  resolve,  and  hard  to  distinguish  from  target  E.  Target  B 
and  C  have  the  same  lateral  position  x  and  y,  and  different  depths.  Target  A  is  close  to  the  source 
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plane,  while  F  is  close  to  the  detector  plane.  As  detailed  in  Section  4  of  Appendix  4,  TROT 
formalism  could  detect  all  6  targets  and  retrieved  their  locations  to  be  at  the  exact  known  target 
locations.  With  the  highly  encouraging  result  from  simulation  even  for  a  considerably  challenging 
task,  we  proceeded  to  test  the  approach  for  the  realistic  situation  of  detecting  and  locating  targets 
from  experimental  data. 

We  have  experimentally  investigated  the  efficacy  of  TROT  by  imaging  both  absorptive  and 
scattering  target(s)  embedded  in  Intralipid-20%  suspension  in  water,  a  model  medium  whose 
optical  absorption  and  scattering  properties  can  be  adjusted  by  varying  the  concentration.  The 
initial  results  for  the  absorptive  target(s)  were  presented  in  the  Third  Annual  Report  [3],  and  are 
detailed  in  Appendix  4.  Here  we  present  a  brief  account  of  the  study  with  a  scattering  target  and 
leave  the  details  for  Appendix  4. 

The  sample  used  a  250  mm  x  250  mm  x  60  mm  transparent  plastic  container  filled  with 
Intralipid-20%  suspension  in  water  as  the  background  medium.  The  concentration  of  Intralipid- 
20%  was  adjusted  to  provide  an  estimated  [8,9]  absorption  coefficient  pa  ~  0.003  mm"1  at  790  nm, 
and  a  transport  mean  free  path  l,  ~  1  mm,  which  were  similar  to  the  average  values  of  those 
parameters  for  human  breast  tissue,  while  the  cell  thickness  of  60  mm  was  comparable  to  thickness 
of  a  typical  compressed  breast.  The  target  was  a  glass  sphere  of  diameter  10  mm  filled  with 
Intralipid-20%  suspension  in  water  to  provide  a  transport  mean  free  path  l,  of  0.25  mm,  and  a 
scattering  coefficient  jus  ~  1 1  mm"1.  The  depth  of  the  scattering  target  was  varied  to  explore  the 
efficacy  of  TROT  in  locating  and  characterizing  a  scattering  target. 

A  multi-source  interrogation  and  multi-detector  signal  acquisition  scheme,  shown  in  Fig.  2  of 
Appendix  4,  was  used  to  acquire  data.  A  100-mW  790-nm  diode  laser  beam  was  used  to  illuminate 
the  samples.  A  1024  x  1024  pixels  charge  coupled  device  (CCD)  camera  equipped  with  a  60-mm 
focal-length  camera  lens  was  used  on  the  opposite  side  of  the  sample  to  detect  the  transmitted  light 
on  the  boundaries  of  the  slab  samples  (detector  plane).  The  pixel  size  was  24  pm.  The  multi-source 
illumination  scheme  was  realized  by  scanning  the  sample  across  the  laser  beam  in  a  two- 
dimensional  x-y  array  of  grid  points  using  a  computer-controlled  translation  stage.  The  sample  was 
scanned  across  the  laser  beam  in  an  array  of  9  x  9  grid  points,  with  a  step  size  of  5  mm.  The 
scanning  and  data  acquisition  processes  were  controlled  by  a  personal  computer  (PC). 

The  depths  (z-positions)  of  the  target  that  were  used  in  the  experiment  are:  15  mm,  20  mm,  25 
mm,  30  mm,  35  mm,  40  mm,  and  45  mm.  A  typical  cross-section  pseudo  image  and  the 
corresponding  spatial  profiles  are  displayed  in  Fig.  1(a)  when  z  =  30  mm.  It  is  compared  to  the 
simulation  results  with  20%  Gaussian  noise  (Fig.  1(b)).  The  lateral  (x,  y)  spatial  profiles  of  the 
pseudo  image  generated  using  simulated  data  are  considerably  wider,  while  the  axial  (z)  spatial 
profile  is  narrower  than  those  obtained  using  experimental  data,  while  the  peak  values  from  the  two 
cases  are  of  the  same  order.  The  retrieved  target  positions  are  listed  in  Table  I. 

It  is  evident  from  the  results  in  Table  I,  that  the  TROT  formalism  locates  the  positions  of  the 
scattering  target  with  considerable  accuracy.  The  lateral  (x,  y)  positions  are  determined  with  higher 
accuracy  than  the  axial  (z)  position.  The  resolution  of  experimental  result  seems  to  be  better  than 
the  simulation,  which  we  attribute  to  the  noise  level  of  20%  used  in  simulation  which  presumably 
is  higher  than  the  experimental  error. 

A  comparison  of  experimental  results  for  scattering  and  absorptive  targets  validate  the  common 
notion  that  it  is  more  challenging  to  locate  and  image  scattering  targets  than  absorptive  targets  in  a 
highly  scattering  medium.  The  results  from  experiment  and  simulation  show  that  TROT,  being  a 
non-iterative  approach,  is  a  faster  and  less  computation  intensive  approach  for  detecting  small 
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targets  in  highly  scattering  turbid  media  and  determining  their  locations  in  3-D  than  other  inverse 
image  reconstruction  techniques. 


Fig.l.  Pseudo  image  of  the  target  (left  pane)  and  corresponding  spatial  intensity  profiles 
(right  pane)  when  the  target  is  located  at  z  =  30  mm:  (a)  experimental  data;  (b)  simulation 
with  20%  Gaussian  noise  added.  P  is  pseudo  value  calculated  using  Eq.  (20)  in  Appendix  4. 


Table  I.  Positions  of  a  scattering  target  located  at  different  depths 


Known 

Positions 
[x,  y,  z  (mm)] 

Retrieved 

Positions 
\x,  y ,  z  (mm)] 

Error 

(mm) 

[Ac,  Ay,  Az\ 

25.7, 24.5,  15 

24.9,  25.9,  18.5 

0.8,  1.4,  3.5 

25.7,  24.5,  20 

27.2,  26.7,  20.5 

1.5,  2.2,  0.5 

25.7, 24.5,  25 

25.7,  26.7, 23.5 

0.0, 2.2,  1.5 

25.7, 24.5,  30 

24.9,  25.2,  32.5 

0.8,  0.7,  2.5 

25.7,  24.5,  35 

24.9,  25.2,  36.5 

0.8,  0.7,  1.5 

25.7,  24.5, 40 

24.9,  25.9,41.5 

0.8,  1.4,  1.5 

25.7, 24.5, 45 

24.9,  25.9, 45.5 

0.8,  1.4,  0.5 

We  are  pursuing  development  of  the  TROT  formalism  for  extended  targets,  which  will  be  suited  to 
detect,  locate,  and  diagnose  breast  tumors  at  later  stages,  as  well. 
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5.  3.  Finite  Element  Method  for  Optical  Tomography 

We  are  exploring  the  Finite  Element  Method  (FEM)  that  has  found  considerable  use  in  optical 
tomography  [13],  as  another  viable  approach  to  optical  mammography  (, Specific  Aim  4,  Task#  17  and  Task# 
18).  One  of  the  objectives  is  to  use  FEM  to  obtain  optical  properties  around  the  suspect  sites  that  TROT  can 
locate  with  high  accuracy  without  needing  long  computation  time.  While  FEM  is  more  computation 
intensive  than  TROT,  using  it  only  over  the  limited  suspect  sites  located  by  TROT  will  reduce  the 
computation  time  significantly. 

We  pursued  testing  and  adaptation  of  a  program  called  NIRFAST  [14]  developed  by  researchers  at 
Dartmouth  College  for  modeling  NIR  frequency  domain  light  transport  in  tissue  based  on  FEM.  We  have 
evaluated  the  program  in  simulation  under  different  conditions  that  include  number  of  targets,  their 
location,  size  and  optical  properties,  sample  geometry,  source  and  detector  positions,  and  noise  level.  It  is 
tested  for  both  two-dimensional  (2 -D)  and  three-dimensional  (3-D)  absorptive  and  scattering  targets,  in 
2D  and  3D  problems.  Some  examples  are  shown  as  below. 

Two-Dimensional  Problems 

First,  we  considered  the  sample  to  be  10-cm  diameter  circle ,  with  background  optical  properties:  fia  = 

0.01  mm"1,  lt=  1  mm  (jus  =  1  mm"1).  Sixteen  (16)  optical  fibers  are  assumed  to  be  placed  at  the  edge 
around  the  circle,  equally  spaced  and  used  as  both  sources  and  detectors.  In  forward  model,  the  source 
positions  are  considered  to  be  located  1  mm  below  the  surface  [15],  where  the  initial  scatterings  of  the 
incident  photons  are  assumed  to  occur.  To  avoid  the  strong  boundary  effect,  when  one  fiber  is  used  as  the 
source,  other  15  fibers  are  used  as  detectors.  Measurements  are  taken  with  every  fiber  used  in  turn  as  a 
source  resulting  in  accumulation  of  data  using  240  source  and  detector  pairs  located  1  mm  below  the  edge 
around  the  circle.  Two  absorptive  and  two  scattering  circular  targets  are  embedded  as  follows. 

Two  absorptive  and  two  scattering  circular  targets  with  different  optical  properties  and  sizes  are 
embedded.  The  center  of  the  circle  is  set  to  be  the  origin  of  the  coordinates  (0,  0)  mm.  The  position, 
optical  property  and  size  of  the  four  targets  are  listed  in  Table  II. 

Table  II.  Target  position,  property  and  size 


Target# 

Position  (x,  y)  (mm) 

Ha  (mm-1) 

Hs'  (mm1) 

Radius  (mm) 

1 

-20,  -20 

0.015 

1 

5 

2 

20,  20 

0.02 

1 

2.5 

3 

20,  -20 

0.01 

2 

5 

3 

-20,  20 

0.01 

4 

2.5 

The  mesh  of  the  sample  and  the  source  (*),  detector  (detector)  and  target  (circle)  positions  are  shown  in 
Fig.  2.  Forward  model  data  in  frequency  domain  with  modulation  frequency  100  MHz  was  generated 
using  NIRFAST.  Background  random  noise  level  of  1%  was  added  and  compared  with  zero-noise  case. 
The  data  is  then  used  for  reconstruction  also  using  NIRFAST  through  an  iterative  process,  where  the 
forward  model  data  is  repeatedly  calculated  while  varying  the  optical  properties  of  every  voxel,  and 
compared  with  the  original  data.  The  uniform  background  property  (jua  =  0.01  mm"1,  jus>  =  1  mm"1)  was 
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Fig.  2.  Sample  mesh  and  the  source, 
detector,  target  positions 
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used  as  the  starting  point.  The  optimization  of 
the  property  is  approached  by  minimizing  the 
difference  between  the  reconstructed  forward 
model  data  and  the  original  data.  Maps  of 
absorption  and  scattering  properties  are 
reconstructed  simultaneously,  as  shown  in  Fig.  3. 

As  shown  in  Fig.  3,  all  four  targets  are  detected, 
and  both  absorption  and  scattering  properties  are 
reconstructed.  The  reconstructed  target  sizes  are  comparable  to  the  actual  ones.  When  no  noise  is  added, 
the  reconstructed  optical  properties  of  the  large  targets  are  close  to  the  known  values,  while  the  properties 
of  the  small  targets  are  not  quite  accurate.  While  1%  random  noise  was  added,  the  sizes  of  the  large 
targets  are  still  close  to  actual  values,  and  the  error  in  reconstructed  properties  is  about  20%.  However,  the 
small  targets  are  hardly  detected.  The  reconstruction  result  can  be  improved  by  changing  conditions,  such 
as,  higher  modulation  frequency,  more  source-detector  pairs  etc. 


Fig.  3.  Reconstructed  absorption  and  scattering  properties 


Next  we  considered  a  rectangular  sample  of  size  80  mm  x  40  mm,  with  /ia  =  0.01  mm'1,  lt  =  1  mm  (jis> 

1  mm'1).  13  sources  and  13  detectors  are  located  on  the  entrance  side  and  exit  side,  respectively.  The 
sources  are  considered  to  be  located  1  mm  below  the  side,  and  used  one  at  a  time  successively. 
Transmitted  data  are  detected  by  13  detectors  for  all  13  sources.  In  total,  169  source  and  detector  pairs  are 
used.  One  absorptive  and  one  scattering  circular  targets  are  embedded.  The  center  of  the  rectangle  is  set  to 
be  the  origin  of  the  coordinates.  The  position,  property  and  size  of  the  targets  are  listed  in  Table  III. 


Table  III.  Target  position,  property  and  size 


Target# 

Position  (x,  y)  (mm) 

Ha  (mm'1) 

HS'  (mm'1) 

Radius  (mm) 

1 

-10,0 

0.02 

1 

5 

2 

10,0 

0.01 

2 

5 
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The  mesh  of  the  sample  with  the  source  (*),  detector  (+)  and  target  (o)  positions  is  shown  in  Fig.  4(a). 
Frequency-domain  data  is  generated  with  modulation  frequency  100  MHz,  and  1%  random  noise  was 
added  and  compared  to  the  situation  with  no  noise. 


NIRFAST  was  then  run  using  the  forward  model  data.  The  reconstructed  images  using  the  absorption  and 
scattering  property  maps  are  shown  in  Fig.  4(b). 


M' 


a 


1%  Noise 


0.009  0.01  0.011  0.012  1  1.2  1.4 

Fig.4(b).  Reconstmcted  absorption  and  scattering  properties 
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As  shown  in  Fig.  4(b),  both  absorptive  and  scattering  targets  were  detected,  and  both  absorption  and 
scattering  properties  were  reconstructed.  The  reconstructed  target  sizes  are  comparable  to  the  actual  sizes 
when  no  noise  is  added.  While  1%  random  noise  is  present,  the  reconstructed  size  is  larger.  When  no 
noise  was  added,  the  error  in  the  reconstructed  optical  properties  is  about  25%  for  the  absorptive  target, 
and  about  15%  for  the  scattering  target.  When  1%  random  noise  was  added,  the  error  in  reconstructed 
properties  is  about  35%  for  the  absorptive  target,  and  about  30%  for  the  scattering  target. 

Three-dimensional  problems 

Next  we  considered  three-dimensional  absorptive  and  scattering  targets.  The  first  sample  was  a  cylinder 
of  diameter  86  mm,  and  height  60  mm,  jua  =  0.01  mm'1,  lt=  1  mm  (jus>=  1  mm'1).  240  source  and  detector 
pairs  were  used,  and  located  in  a  manner  similar  to  the  2-D  circular  problem  discussed  in  the  previous 
section  in  the  cross-section  of  the  middle  plane  (z  =  0  mm)  of  the  cylinder.  One  absorptive  (scattering) 
spherical  target  was  embedded  at  (x,  y,  z)  =  (20,  0,  0)  mm.  A  discretized  mesh  of  the  cylindrical  sample  is 
shown  in  Fig.  5,  along  with  the  positions  of  the  target  (dark  red),  sources  (green  circle)  and  detectors 
(blue  plus).  Frequency-domain  forward  model  data  was  generated  with  modulation  frequency  100MHz 
and  1%  random  noise  added  compared  with  no  noise  present  using  NIRFAST,  and  then  fed  back  into  to 
reconstruct  the  absorption  and  scattering  property  maps  of  the  sample. 


£ 

E 


Fig. 5.  Sample  mesh  (left  pane)  and  source/detector/target  positions  (right  pane) 


Absorptive  target 

An  absorptive  target  is  embedded,  with  =  0.1  mm'1,  jus  the  same  as  background,  and  radius  5  mm. 
Reconstructed  3D  image  and  a  cross-section  through  the  middle  of  the  target  are  shown  in  Fig.  6(a)  when 
there  is  no  noise,  and  in  Fig.  6(b)  with  1%  noise. 
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Figure  6(a).  Reconstructed  images  (absorption  and  scattering  properties)  of  an  absorptive  target  under  no 
background  noise  condition:  (top)  3-D  distribution  of  absorption  coefficient;  (bottom  left)  absorption  coefficient 
distribution,  and  (bottom  right)  reduced  scattering  coefficient  distribution  through  the  cylinder  center  (z  =  0)  plane. 
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1%  noise: 
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Figure  6(b).  Reconstructed  images  (absorption  and  scattering  properties)  of  an  absorptive  target  with  1% 
background  noise:  (top)  3-D  distribution  of  absorption  coefficient;  (bottom  left)  absorption  coefficient  distribution, 
and  (bottom  right)  reduced  scattering  coefficient  distribution  through  the  cylinder  center  (z  =  0)  plane. 


Both  absorptive  and  scattering  properties  are  reconstructed  simultaneously.  Even  though  it  is  an 
absorptive  target,  and  the  scattering  property  of  the  target  is  the  same  as  that  of  the  background,  the  target 
is  also  seen  in  the  reconstructed  scattering  property  map.  The  reconstructed  target  size  is  comparable  to 
the  actual  size  in  the  lateral  direction,  while  much  bigger  in  the  z  direction  for  both  no  noise  and  1%  noise 
level  cases.  The  reconstructed  absorptive  property  of  the  target  is  about  15%  of  the  actual  value.  The 
results  can  be  improved  by  using  higher  modulation  frequency,  and  more  source-detector  pairs  etc. 
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Scattering  target 


A  scattering  target  of  radius  5  mm  with  jus  =  10  mm"1,  and  jua  the  same  as  background  is  embedded  and 
radius  5  mm.  The  3D  absorption  and  scattering  characteristics  maps  were  reconstructed  in  a  manner 
similarly  to  that  used  for  the  absorptive  target.  The  reconstructed  3D  image  and  a  cross-section  through 
the  middle  of  the  target  are  shown  in  Fig.  7(a)  when  the  background  noise  level  is  considered  to  be  0  and 
in  Fig.  7(b)  when  1%  noise  is  assumed  to  be  present. 

No  noise 

Ps 


y(mm)  y(mm) 

Fig.  7(a).  Reconstructed  scattering  and  absorption  properties  of  the  single  scattering  target  assuming  0  background 
noise  (top)  3-D  distribution  of  reduced  scattering  coefficient;  (bottom  left)  absorption  coefficient  distribution,  and 
(bottom  right)  reduced  scattering  coefficient  distribution  through  the  cylinder  center  (z  =  0)  plane. 
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Fig  7(b)  Reconstructed  scattering  and  absorption  properties  of  the  single  scattering  target  assuming  1%  background 
noise:  (top)  3-D  distribution  of  reduced  scattering  coefficient;  (bottom  left)  absorption  coefficient  distribution,  and 
(bottom  right)  reduced  scattering  coefficient  distribution  through  the  cylinder  center  (z  =  0)  plane. 

Both  absorptive  and  scattering  properties  were  reconstructed  simultaneously  in  this  simulation.  Even 
though  it  is  a  scattering  target,  it  appears  in  the  reconstructed  absorption  property  map  as  well,  especially 
for  no-noise  case.  The  reconstructed  target  size  is  comparable  to  the  actual  size  in  the  lateral  direction  for 
both  no  noise  and  1%  noise  level  cases.  The  reconstructed  scattering  property  of  the  target  is  about  11- 
12%  of  the  actual  value  for  no  noise  and  1%  noise  level  cases. 

Two  targets 

The  next  situation  involved  two  absorptive  spherical  targets  embedded  at  (-20,  0,  0)  mm  and  (20,  0,  0) 
mm  of  the  sample  volume.  Each  of  the  spheres  has  a  radius  of  5  mm,  fia  =  0.1  mm"1  and  jus  is  the  same  as 
background.  The  cylinder  mesh  and  the  source/detector  positions  and  background  optical  characteristics 
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of  the  sample  are  the  same  as  in  the  single  target  case  reported  above.  Frequency-domain  data  were 
generated  with  modulation  frequency  of  100  MHz.  Figure  8(a)  shows  the  reconstructed  absorption  and 
scattering  coefficient  maps  of  the  target  for  0  background  noise,  while  Figure  8(b)  shows  those  when  1% 
background  noise  is  assumed.  A  cross-section  image  through  the  center  of  the  cylinder  is  also  shown. 

No  noise: 

Ma 


Fig.  8(a)  Reconstmcted  scattering  (and  absorption)  properties  of  the  targets  assuming  0  background  noise:  (top)  3-D 
distribution  of  absorption  coefficient;  (bottom  left)  absorption  coefficient  distribution,  and  (bottom  right)  scattering 
coefficient  distribution  through  the  cylinder  center  (z  =  0)  plane. 
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1%  noise: 
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Fig.  8(b).  Reconstructed  scattering  (and  absorption)  properties  of  the  targets  assuming  1%  background  noise: 

(top)  3-D  distribution  of  absorption  coefficient;  (bottom  left)  absorption  coefficient  distribution,  and  (bottom  right) 
scattering  coefficient  distribution  through  the  cylinder  center  (z  =  0)  plane. 


Overall,  the  target  size  can  be  reconstructed  with  values  comparable  to  the  actual  value  for  both  2D  and 
3D  cases.  It  is  much  more  difficult  to  reconstruct  optical  property  of  the  target(s)  accurately.  What  is  even 
more  noteworthy  is  that  even  for  an  absorptive  (scattering)  target  a  reduced  scattering  (absorption) 
coefficient  map  was  obtained  in  addition  to  the  desired  absorption  (reduced  scattering)  coefficient  map. 
This  is  not  desirable  and  needs  to  be  addressed.  Higher  modulation  frequency  and  more  source-detector 
pairs  may  be  used  to  improve  the  results. 

The  program  is  now  being  tested  for  slab  geometry  which  is  used  in  our  project.  We  plan  to  pursue  the 
FEM  approach  further  for  3D  image  reconstruction  and  estimation  of  target  optical  properties. 
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5.4  Magnetic  Resonance  Spectroscopic  Imaging  Training 

The  CCNY  researchers  pursued  receiving  training  ( Specific  Aim  0,  Tasks  2,  3  and  5)  on  magnetic 
resonance  spectroscopic  imaging  (MRSI)  in  primary  collaborating  mentor  (PCM)  Dr.  Koutcher’s 
lab  at  MSKCC.  The  thrust  of  the  effort  was  to  adapt  MRSI  approach  for  detection  of  lactate, 
which  is  expected  to  provide  a  window  to  explore  tumor  aggressiveness  (Specific  Aim  1).  Once 
MRSI  is  established  as  a  reference  technique,  the  aim  is  to  develop  optical  spectroscopic 
techniques  for  lactate  detection  and  validate  the  measurements  against  MRSI  results.  Since 
lactate  level  is  associated  with  tumor  aggressiveness,  metastasis  and  treatment  response  [16], 
dependable  and  noninvasive  means  for  detecting  lactate  in  tumors  is  of  immense  interest. 

Conventional  MRSI  is  not  adequate  for  effective  lactate  detection  since  the  overlapping  lipid  and 
adjacent  water  signals  obscure  lactate  signal.  The  double  frequency-selective  multiple  quantum 
coherence  transfer  (SelMQC)  technique  shows  promise  to  detect  lactate  [17],  and  the  trainees  got 
involved  in  developing  instrumentation  and  testing  those  for  lactate  detection  using  phantom  in 
Year  3  of  the  project.  They  pursued  their  training  further  during  this  reporting  period.  We  plan  to 
undertake  correlated  MRSI  with  SelMQC  and  optical  spectroscopic  measurements  to  investigate 
the  glycolytic  activity  of  tumors,  and  explore  their  potential  for  tumor  aggressiveness 
assessment. 

5.5  Research  Proposal  Development 

One  of  the  proposed  tasks  ( Specific  Aim  0,  Task#  7)  involves  developing  at  least  one  research 
proposal  and  submitting  it  to  NIH  or  USAMRMC  for  funding.  We  submitted  a  pre-proposal 
entitled,  “ Multi-functional  tumor-targeting  nanocomposites  and  time-reversal  optical  imaging 
for  early  detection  of  breast  cancer  and  prevention  of  micro-metastases ”  to  the  Idea  Award 
(Collaborative  Option)  category  of  the  2011  Breast  Cancer  Research  Program  of  CDMRP.  S.  K. 
Gayen  (PI  of  this  proposal)  is  the  Initiating  PI  of  the  above-mentioned  proposal.  The 
Collaborative  PI  is  Valeria  Balogh-Nair  of  the  CCNY  Chemistry  Department.  A  brief  overview 
of  the  pre-proposal  follows: 

Research  Idea 

The  objective  of  the  proposed  research  is  to  develop  multi-functional,  tumor- targeting 
nanocomposites  and  near-infrared  (NIR)  optical  imaging  approaches  for  early  detection  of  breast 
cancer,  and  prevention  of  micro-metastases  that  are  responsible  for  majority  of  breast  cancer 
mortality.  The  nanocomposite  synthesis  will  use  a  multivalent  dendrimer  (a  class  of  organic 
macromolecules)  platform  to  incorporate  fluorescent  moieties  (e.g.,  semiconductor  quantum 
dots,  or  gold/silver  nanoparticles)  as  contrast  agents  for  imaging,  and  chemokine  mimics  for 
selective  targeting  of  cancer  cells  and  prevention  of  metastases.  The  high  affinity  of  chemokine 
mimic’s  ligands  for  chemokine  CXCR4  receptors  will  ensure  selective  delivery  of  the 
nanocomposite  to  the  target.  A  multi-source  NIR  probing  and  multi-detector  signal  acquisition 
arrangement  along  with  a  time-reversal  image  reconstruction  algorithm  will  be  used  for  fast 
detection  of  tumor  and  determination  of  its  location  in  three  dimensions. 
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Impact  and  Innovation 

The  major  impact  of  the  proposed  research  is  that  it  has  the  potential  to  provide  a  modality 
for  early  detection  of  breast  cancer,  and  for  prevention  of  metastases,  the  major  cause  of 
mortality.  A  broader  potential  impact  is  that  the  approach  could  be  extended  to  other  cancers 
using  chemokine  mimics  directed  to  other  chemokine  receptors. 

The  project  is  highly  innovative  in  many  ways.  First,  the  synthesis  process  brings  together 
the  concept  of  multivalency,  the  idea  of  chemokine  mimics  for  prevention  of  cancer  cell 
migration,  and  use  of  dendrimers  as  stabilizers  and  nanoreactors  for  contrast  agent  synthesis.  A 
combination  of  these  novel  ideas  is  proposed,  to  the  best  of  our  knowledge  for  the  first  time,  as  a 
multi-prong  approach  to  fight  the  menace  of  breast  cancer.  Second,  the  affinity  of  the  chemokine 
mimics  to  seven  helix  chemokine  CXCR4  receptors  obviates  the  need  for  a  targeting  vector.  The 
chemokine  mimics  will  play  the  dual  role  of:  (1)  “homing  devices”  for  selective  delivery  of  the 
nanocomposite  to  the  tumor  sites;  and  (2)  “prevention  agents”  interacting  with  the  chemokine 
receptor  sites  to  inhibit  metastases.  Third,  the  efficacy  of  noninvasive  NIR  imaging  approach  for 
early  detection  and  potential  diagnosis  will  be  significantly  enhanced  through  design  of  efficient 
contrast  agent.  Fourth,  the  idea  of  time -reversal  optical  tomography  (TROT)  is  a  new  paradigm 
in  diffusive  optical  tomographic  (DOT)  imaging.  While  currently  pursued  DOT  approaches  are 
iterative  and  computation  time  intensive,  TROT  is  non-iterative  and  will  be  faster,  which  is  a 
necessary  condition  for  real-time  imaging.  TROT  is  designed  for  detecting  and  locating  small 
targets  and  will  be  suited  for  early  detection  when  tumors  are  small.  Finally,  the  dendrimer 
platform  with  multivalent  surface  will  enable  incorporation  of  moieties  that  enhance  other 
imaging  modalities,  such  as,  magnetic  resonance  imaging,  enabling  development  of  sought-after 
dual  or  multimodal  imaging  modules. 

6.  KEY  ACCOMPLISHMENTS 

•  Key  research  accomplishments  include:  (a)  development  of  Time  Reversal  Optical 
Tomography  (TROT)  (which  is  fast  and  minimally  iterative),  and  demonstration  of  its 
efficacy  in  detecting  small  targets  in  a  tissue  simulating  turbid  medium;  (b)  adaptation  of 
decomposition  methods,  such  as,  Principal  Component  Analysis  (PCA)  and  nonnegative 
matrix  factorization  (NMF)  for  study  of  the  diffuse  optical  imaging  problem  and  comparing 
those  with  previously  developed  OPTICA  (optical  tomography  using  Independent 
Component  Analysis);  (c)  exploring  lactate  detection  using  magnetic  resonance 
spectroscopic  imaging  as  a  potential  method  for  assessing  tumor  aggressiveness;  (d) 
publication  of  several  papers  and  presentation  of  research  results  in  major  conferences 
including  Era  of  Hope  (2011)  (detailed  in  Section  7  below);  and  (e)  investigating  finite 
element  method  with  the  aim  of  using  it  as  a  complement  to  TROT  for  retrieval  of  optical 
parameters  at  the  localized  target  sites. 

•  The  key  training  accomplishment  includes  the  successful  participation  of  physical  scientists 
and  engineers  of  CCNY  research  team  in  cancer  biology  research  involving  magnetic 
resonance  spectroscopic  imaging  at  the  MRI  research  facility  of  MSKCC. 

•  Development  and  submission  of  an  independent  research  pre-proposal  to  BCRP  2011  is 
indicative  of  the  CCNY  team’s  progress  towards  developing  a  sustainable  breast  cancer 
research  program  at  CCNY. 
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7.  REPORTABLE  OUTCOMES 
Publications 

(1)  B.  Wu,  W.  Cai,  M.  Alrubaiee,  M.  Xu  and  S.  K.  Gayen,  “ Three-dimensional  time-reversal 
optical  tomography in  Multimodal  Biomedical  Imaging  VI.  Proceeding  of  SPIE,  vol. 
7892,  pp.  78920G-1  -  78920G-6  (2011).  ( Attached  as  Appendix  1) 

(2)  M.  Alrubaiee,  B.  Wu,  M.  Xu,  W.  Cai,  and  S.  K.  Gayen,  “Multi-wavelength  diffusive 
optical  tomography  using  Independent  Component  Analysis  and  Time  Reversal 
algorithms”,  in  Diffuse  Optical  Imaging  III  Proceeding  of  SPIE-OSA,  vol.  8088,  pp. 
80880Y-1  -  80880Y-5  (2011).  ( Attached  as  Appendix  2) 

(3)  Binlin  Wu,  M.  Alrubaiee,  W.  Cai,  M.  Xu  and  S.  K.  Gayen,  “Diffuse  Optical  Imaging 
using  Decomposition  Methods.”  International  Journal  of  Optics  (Accepted  for 
publication)  {Attached  as  Appendix  3) 

(4)  Binlin  Wu,  W.  Cai,  M.  Alrubaiee,  M.  Xu,  and  S.  K.  Gayen,  “Time  reversal  optical 
tomography:  locating  targets  in  a  highly  scattering  turbid  medium,”  Opt.  Express  19, 
21956-21976  (2011).  {Attached  as  Appendix  4) 

Conference  Presentations 

(1)  M.  Alrubaiee,  Binlin  Wu,  W.  Cai,  M.  Xu  and  S.  K.  Gayen,  “Multi-wavelength  diffusive 
optical  tomography  using  independent  component  analysis  and  time  reversal  algorithms.” 
Paper  8088-33  presented  at  the  SPIE/OSA  European  Conference  on  Biomedical  Optics, 
22-26  May,  2011,  Munich,  Germany. 

(2)  Binlin  Wu,  W.  Cai,  M.  Alrubaiee,  M.  Xu  and  S.  K.  Gayen,  "Three  dimensional  time 
reversal  optical  tomography,"  Paper  7892-15  presented  at  the  SPIE's  International 
Symposium  on  Biomedical  Optics,  BiOS  '11/  Photonics  West,  22-27  January,  San 
Francisco,  California. 

(3)  Binlin  Wu,  W.  Cai,  M.  Alrubaiee,  M.  Xu,  and  S.  K.  Gayen,  “Time  reversal  optical 
tomography  for  breast  tumor  detection”,  Paper  BC060 114-2962  presented  at  the  6th  Era 
of  Hope  2011  conference  of  Department  of  Defense  (DOD)  Breast  Cancer  Research 
Program  (BCRP),  2-5  August,  Orlando,  Florida. 

Research  Proposal 

(1)  S.  K.  Gayen  (Initiating  PI),  V.  Balogh-Nair  (Collaborating  PI),  “ Multi-functional  tumor¬ 
targeting  nanocomposites  and  time-reversal  optical  imaging  for  early  detection  of  breast 
cancer  and  prevention  of  micro-metastases submitted  to  the  Idea  Award  (Collaborative 
Option)  category  of  the  201 1  Breast  Cancer  Research  Program  of  CDMRP. 

Research  Personnel  Development 

Xiaohui  Ni,  a  researcher  trained  and  supported  in  part  by  this  grant  has  moved  to  the 

Department  of  Chemistry  and  Chemical  Biology  of  Harvard  University  as  a  research 

associate. 

8.  CONCLUSION 

The  work  carried  out  during  this  reporting  period:  (a)  culminated  in  CCNY  research  trainees’ 
participation  in  magnetic  resonance  spectroscopic  imaging  breast  cancer  research  at  MSKCC; 
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and  (b)  shows  the  potential  for  noninvasive  detection  and  three-dimensional  localization  of  a 
tumor  within  a  breast  with  significant  accuracy.  The  contrast  is  based  on  the  differences  in  the 
light  scattering  and  absorption  characteristics  of  the  tumor  and  normal  breast  tissue. 

“So  What  Section” 

•  The  National  Cancer  Institute  (NCI)  has  identified  the  development  of  imaging 
methodologies  as  an  extraordinary  opportunity  for  advancement  in  cancer  research.  Since  the 
background  of  the  CCNY  team  is  in  physical  sciences  and  engineering,  the  training  they 
received  has  provided  them  with  necessary  laboratory  background  in  the  biology  of  cancer 
research,  and  helping  develop  a  knowledgeable  multidisciplinary  research  force  in  the  fight 
against  breast  cancer. 

•  A  recent  study  involving  35,319  patients  underscores  the  influence  of  primary  tumor  location 
on  breast  cancer  prognosis  [18],  and  makes  it  imperative  that  breast  cancer  detection 
modalities  to  obtain  three-dimensional  (3-D)  location  of  the  tumor  relative  to  the  axilla  be 
developed.  The  optical  imaging  techniques  (TROT,  OPTICA  and  other  decomposition 
methods)  being  developed  are  an  important  steps  in  obtaining  3-D  location  of  a  tumor  within 
the  breast.  The  methods  are  minimally  iterative,  fast,  and  designed  for  locating  small  targets, 
which  make  those  suitable  for  detecting  tumors  in  early  stages  of  development  when  those 
are  more  amenable  to  treatment. 

•  The  study  involving  lactate  detection  holds  promise  for  tumor  aggressiveness  assessment. 
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ABSTRACT 

Time  reversal  optical  tomography  (TROT)  approach  is  used  to  detect  and  locate  absorptive  targets  embedded  in  a  highly 
scattering  turbid  medium  to  assess  its  potential  in  breast  cancer  detection.  TROT  experimental  arrangement  uses  multi¬ 
source  probing  and  multi-detector  signal  acquisition  and  Multiple- Signal-Classification  (MUSIC)  algorithm  for  target 
location  retrieval.  Light  transport  from  multiple  sources  through  the  intervening  medium  with  embedded  targets  to  the 
detectors  is  represented  by  a  response  matrix  constructed  using  experimental  data.  A  TR  matrix  is  formed  by  multiplying 
the  response  matrix  by  its  transpose.  The  eigenvectors  with  leading  non-zero  eigenvalues  of  the  TR  matrix  correspond  to 
embedded  objects. 

The  approach  was  used  to:  (a)  obtain  the  location  and  spatial  resolution  of  an  absorptive  target  as  a  function  of  its 
axial  position  between  the  source  and  detector  planes;  and  (b)  study  variation  in  spatial  resolution  of  two  targets  at  the 
same  axial  position  but  different  lateral  positions.  The  target(s)  were  glass  sphere(s)  of  diameter  ~9  mm  filled  with  ink 
(absorber)  embedded  in  a  60  mm-thick  slab  of  Intralipid-20%  suspension  in  water  with  an  absorption  coefficient  fia  ~ 
0.003  mm"1  and  a  transport  mean  free  path  ~  1  mm  at  790  nm,  which  emulate  the  average  values  of  those  parameters 
for  human  breast  tissue.  The  spatial  resolution  and  accuracy  of  target  location  depended  on  axial  position,  and  target 
contrast  relative  to  the  background.  Both  the  targets  could  be  resolved  and  located  even  when  they  were  only  4-mm 
apart.  The  TROT  approach  is  fast,  accurate,  and  has  the  potential  to  be  useful  in  breast  cancer  detection  and  localization. 

Keywords:  Time  reversal,  MUSIC,  diffusive  optical  imaging,  optical  tomography,  biomedical  imaging,  breast  cancer 
imaging,  scattering  medium,  diffusion  approximation 


1.  INTRODUCTION 

Optical  imaging  of  targets  embedded  in  turbid  media,  such  as  a  tumor  in  a  breast,  has  attracted  much  attention  in  the  last 
two  decades.  When  a  beam  of  light  propagates  through  a  highly  scattering  medium,  photons  are  scattered  and  diffused 
into  a  broad  area;  phase  coherence  and  polarization  of  light  deteriorate;  short  pulses  broaden;  and  consequently  sharp 
images  of  the  targets  cannot  be  formed  directly.  Various  algorithms  have  been  developed  to  perform  image 
reconstruction1'4.  Inverse  image  reconstruction  (HR)  is  an  ill-posed  problem  and  search  of  reliable  and  fast  approaches  is 
an  important  and  formidable  task  for  optical  imaging  of  human  tissue  such  as  breast.  Recent  inverse  algorithms,  such  as 
Newton-Raphson-Marquart  algorithms5  and  direct  linear  inversion  of  3-D  matrices6,  are  time  consuming.  The  iterative 
methods6, 7  may  not  ensure  that  the  obtained  result  arrives  at  a  “global  minimum”  or  converge  to  a  “local  minimum”. 

Time  reversal  optical  tomography  (TROT)9,10  was  introduced  as  a  reconstruction  method  for  imaging  targets  in 
turbid  media  non-invasively  using  light,  following  the  development  of  TR  imaging  using  multistatic  radar  signal8. 
Compared  to  other  inverse  methods  which  usually  require  iterations4,  TROT  is  fast  since  there  is  no  iteration  involved. 
The  signals  and  noise  are  separated  into  orthogonal  subspaces.  A  method  similar  to  L-curve  regularization  is  used  to 
select  the  signal  subspace.  Then  a  pseudo  spectrum  is  calculated  directly  for  all  voxels  in  the  sample  using  the  vector 
subspace  method,  Multiple  Signal  Classification  (MUSIC)8"10.  Tomographic  pseudo  images  can  be  generated  using 
pseudo  values.  Locations  and  characteristics  of  targets  are  determined  by  the  global  maximum  (or  local  maximum  in  low 
signal-to-noise  (SNR)  cases)  components  in  the  pseudo  spectrum. 
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This  paper  is  organized  as  follows.  In  section  2,  the  formalism  of  the  TROT  approach  is  presented.  In  section  3,  the 
experimental  arrangement,  materials  and  parameters  are  described.  In  section  4,  the  TROT  analysis  on  the  experimental 
data  and  results  are  presented.  Section  5  serves  as  discussion  and  summary. 


2.  FORMALISM 

Using  diffusion  approximation  of  the  radiative  transfer  equation  (RTE),  the  scattered  light  due  to  small  weak 
inhomogeneities  (targets)  embedded  in  a  homogeneous  medium  to  the  first  order  Bom  approximation  can  be  written  as4 

4>sca(rd,rs )  =  -  f  G(rdlr)S[ia(r)cG(r,rs)d3r  -  f  5D(r)cVrG(rd,r)  ■  VrG(r,rs)d3r  ,  (1) 

where  rs,  rd,  and  r  are  the  positions  of  a  point-like  source  of  unit  power,  detector  and  target,  respectively;  G(r,  rs)  and 
G(rd ,  r)  are  the  Green  functions  that  describe  light  propagations  from  the  source  to  the  target  and  from  the  target  to  the 
detector,  respectively;  Sjua  and  3D  describe  the  differences  of  absorption  and  scattering  properties  between  the  targets 
and  the  background  medium,  respectively;  and  c  is  the  light  speed  in  the  medium.  For  an  absorptive  target,  a  matrix  form 
of  the  response  matrix  K  is  constmcted  as  {K^}  =  ( i  =  1,2 ,-mmNd;  j  =  1,2,  ••• , Ns), 

where  rh  Tj  and  Xm  are  locations  of  the  ith  detector,  jth  source  and  mth  target,  respectively;  Ns,  Nd  and  M  are  the 
numbers  of  sources,  detectors  and  targets,  respectively.  It  is  assumed  the  number  of  targets  is  less  than  the  number  of 
sources  and  detectors,  M<  min(Nd ,  Ns). 

A  multi-source  interrogation  and  multi- detector  acquisition  scheme  is  used  to  acquire  multistatic  trans illumination 
data,  from  which  the  scattered  light  due  to  the  targets  is  found,  (j)sca  =  0  —  0O,  where  0  is  the  light  intensity  measured 
on  the  boundary  with  targets  embedded  in  the  scattering  medium  and  0O  is  the  background  image  without  targets 
embedded,  which  can  be  approximated  by  an  “average”  of  all  acquired  images.  Thus,  the  response  matrix  K  is 
constructed  with  perturbations  of  spatial  intensity  distributions  on  the  boundary.11 

A  time  reversal  matrix  T  is  then  constructed  as  T  =  KK ^  (or  T  =  KKT  for  continuous  wave  (CW)),  which  is  similar 
to  the  time  reversal  matrix  used  in  the  time  reversal  by  Devaney8"11.  Eigenvalues  Aj  and  eigenvectors  Uj  of  T  are  found, 
where  j  =  1,2,  •••  min(Nd,  Ns).  The  leading  eigenvalues  correspond  to  the  targets.  The  vector  subspace  method,  MUSIC 
is  used  to  determine  the  locations  of  hidden  objects.  A  Green’s  function  vector  on  the  detector  array  g^Xp)  associated 
with  a  test  target  position  Xp  at  pth  voxel  is  calculated  using  diffusion  model,  where 
#d(Xp)  —  [Gd(ri>Xp),Gd(r2,Xp),---  ,Gd(rNd,Xp)]T  and  the  Green’s  functions  Gd(rj,Xp ),  j  =  1,2 ,-mmNd,  describe 
light  propagation  from  the  test  target  position  Xp  to  detectors  at  r;-.  Then  the  eigenvectors  Uj  corresponding  to  the 
leading  eigenvalues  are  used  to  calculate  the  MUSIC  type  pseudo  spectrum  using  the  following  formula8,9: 

pfy  1  _  _ \3d(Xp)\2 _  m 

{  V)  \\3d(Xp)\2-Z^e\uS9d(Xp)\2\’  ^  ^ 

where  e  is  the  threshold  determined  using  an  L-curve  method  to  separate  the  signal  and  noise  subspaces.  When  both 
absorption  and  scattering  properties  are  considered12,  one  Green’s  function  vector  associated  with  absorption  property 
and  constructed  with  Green’s  functions  Gd  and  three  Green’s  function  vectors  associated  with  scattering  property 
constructed  with  dGd/dx,  dGd/dy  and  dGd/dz ,  respectively  are  used  to  calculate  the  pseudo  spectrum  at  each  voxel. 
Hence,  four  pseudo  values  are  obtained  at  each  voxel,  one  for  absorption,  and  three  for  scattering. 

A  maximum  value  of  E(2fp)  is  obtained  when  Xp  is  the  position  of  one  of  the  hidden  objects.  By  sorting  pseudo 
values  P(Xp),  the  positions  of  the  embedded  objects  are  determined.  At  the  same  time,  absorptive  and  scattering  objects 
are  distinguished  according  to  the  type  of  the  Green’s  function  vector  g^Xp)  associated  with  maximum  pseudo  value. 

In  this  MUSIC  procedure,  only  a  single  run  for  calculating  the  pseudo  spectrum  over  voxels  is  required,  without 
iterative  procedure  used  in  the  traditional  inverse  approach. 
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3.  EXPERIMENTAL  METHODS  AND  MATERIALS 


Two  different  samples  were  used  to  test  the  efficacy  of  the  TROT  approach  in  two  different  experiments.  Both  samples 
used  a  250  mm  x  250  mm  x  60  mm  transparent  plastic  container  filled  with  Intralipid-20%  suspension  in  water.  The 
concentration  of  Intralipid-20%  was  adjusted  to  provide  an  absorption  coefficient  / ia  ~  0.003  mm"1  at  790  nm,  and  a 
transport  mean  free  path  lt~  1  mm,  which  were  similar  to  values  of  those  parameters  for  human  breast  tissue.  In  the  first 
sample,  the  depth  (position  along  z-axis)  of  an  absorptive  target  was  varied  to  explore  how  the  accuracy  of  position 
estimate  depended  on  depth.  In  the  second  sample,  the  separation  between  two  absorptive  targets  was  varied  to  test  how 
close  those  could  be  and  yet  be  resolved  as  separate  objects.  The  target(s)  were  glass  sphere(s)  of  diameter  8-9  mm  filled 
with  ink  dissolved  in  Intralipid  suspension  in  water  (jus  was  adjusted  to  be  the  same  as  that  of  the  background  medium,  fia 
=0.013  mm'1  which  was  about  3  times  higher  than  that  of  background  medium). 

A  multi-source  interrogation  and  multi- detector  acquisition  scheme  was  used  to  acquire  data.  A  100-mW  790-nm 
diode  laser  beam  was  used  to  illuminate  the  samples  and  scan  the  source  plane  with  a  5  mm  step-size.  A  charge  coupled 
device  (CCD)  camera  was  used  on  the  other  side  of  the  sample  to  detect  the  transmitted  light  on  the  boundaries  of  the 
slab  samples  (detector  plane).  The  first  sample  was  scanned  with  the  laser  beam  in  an  array  of  9x9  grid  points,  and  the 
second  sample  was  scanned  in  an  array  of  11x15  grid  points.  The  scanning  and  acquisition  process  was  controlled  by  a 
computer  (PC).  A  schematic  diagram  of  the  experimental  setup  is  shown  in  Fig.  1. 


◄ - PC 


Fig.  1.  A  schematic  diagram  of  the  experimental  arrangement  used  for  imaging  objects  embedded  in  a  turbid  medium. 


4.  ANALYSIS  AND  RESULTS 

One  whole  dataset  for  each  experiment  is  of  order  108  source-detector  pairs,  considering  each  pixel  in  a  CCD  camera  as  a 
detector.  From  each  image,  a  region  of  interest  was  cropped  out  and  then  every  5x5  pixels  in  the  cropped  image  were 
binned  to  one  pixel.  Light  intensity  perturbation  due  to  targets  in  each  image  was  found  by  subtracting  the  background 
image  from  each  individual  image.  The  response  matrix  was  constructed  using  the  light  intensity  perturbations.  The  TR 
matrix  was  generated  by  multiplying  the  response  matrix  by  its  adjoint  matrix  (transpose  for  continuous-wave  (CW) 
illumination).  The  eigenvalue  equation  was  solved  and  signal  subspace  was  selected  and  separated  from  the  noise  subspace. 
MUSIC  was  then  used  to  calculate  the  pseudo  spectrum  for  every  voxel  in  the  3D  space  of  the  sample.  For  each  voxel,  one 
absorption  and  three  scattering  components  were  calculated.  The  voxel  size  was  0.77  mm  x  0.77  mm  x  1  mm.  By  sorting 
the  pseudo  spectrum,  the  object(s)  were  located. 

In  the  first  experiment  using  only  one  target,  the  lateral  (x,  y)  position  of  the  target  was  kept  the  same  at  (25.5  mm,  24.7 
mm),  while  five  different  depths  (position  along  z-axis)  of  15  mm,  20  mm,  25  mm,  30  mm,  35  mm,  40  mm  and  45  mm 
were  used.  A  cross-sectional  pseudo  image  was  generated  using  the  pseudo  spectrum  for  all  voxel  positions  in  the  sample. 
The  spatial  profiles  in  the  x,  y  and  z  directions  of  the  images  through  the  targets  (a  typical  image  and  profiles  for  the  two- 
target  case  are  shown  in  Fig.  2)  are  used  to  assess  target  location.  The  retrieved  target  positions  are  compared  with  known 
positions  in  Table  1 . 


Proc.  of  SPIE  Vol.  7892  78920G-3 


Downloaded  from  SPIE  Digital  Library  on  11  Feb  2011  to  68.173.227.221.  Terms  of  Use:  http://spiedl.org/terms 


Table  1.  Positions  of  one  target  located  at  different  depths 


Known  Positions 

Retrieved  Positions 

Error  (mm) 

[x,y,z(  mm)] 

[x,  y,  z  (mm)] 

25.5,  24.7,  15 

24.9, 24.4, 17.5 

0.6,  0.3,  2.5 

25.5,  24.7,  20 

25.7,  24.4,21.5 

0.2,  0.3,  1.5 

25.5,  24.7,  25 

25.7,  24.4,  26.5 

0.2,  0.3,  1.5 

25.5,  24.7,  30 

25.7,  24.4,  30.5 

0.2,  0.3,  0.5 

25.5,  24.7,  35 

25.7,  25.2,  33.5 

0.2,  0.5,  1.5 

25.5,  24.7,  40 

24.9,  25.2,  36.5 

0.6,  0.5,  3.5 

25.5,  24.7,  45 

24.9,  25.2,  39.5 

0.6,  0.5,  5.5 

The  TROT  assessed  positions  were  in  good  agreement  with  the  known  positions.  The  accuracy  of  the  z-position  was 
found  to  be  optimal  when  the  target  was  located  in  the  middle  plane  of  the  sample,  and  deteriorated  when  the  target  was 
closer  to  the  source  plane  or  the  image  plane. 


In  the  second  experiment  using  two  targets  the  depth  was  kept  fixed  (z  =  30  mm),  while  the  separation  between  them 
separated  with  distances  of  ~  11.8  mm,  16.8  mm,  and  26.8  mm,  in  the  x  direction.  A  cross-sectional  pseudo  image  of  the 
targets  when  separated  by  a  center-to-center  distance  of  11.8  mm  (separation  between  nearest  edges  ~4  mm),  generated 
using  the  pseudo  spectrum  are  shown  in  the  left  pane  of  Fig.  2.  Similar  images  for  other  target  separations  were  also 
obtained.  The  profiles  in  the  x,  y  and  z  directions  through  the  right  target  are  shown  in  the  right  pane  of  Fig.  2.  These 
profiles  were  used  to  assess  locations  of  the  targets,  and  the  separation  between  the  two  targets.  In  all  cases,  the  targets  were 
determined  to  be  absorptive,  because  peaks  occurred  in  the  pseudo  spectrum  with  the  test  Green’s  function  vector 
corresponding  to  absorption  property. 
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Fig.  2.  TROT  generated  cross-section  pseudo  image  when  the  targets  are  separated  by  11.8  mm  is  shown  in  the  left  pane  and  pseudo¬ 
value  profiles  through  the  target  along  x,  y  and  z  directions  are  shown  in  the  right  pane. 


The  known  and  retrieved  positions  and  separations  Ax  between  of  the  two  targets  appear  in  Table  2.  In  all  the  cases,  the 
two  targets  were  resolved,  even  when  their  center-to-center  separation  was  1 1.8  mm  apart,  nearest  sides  separated  by  only 
~4  mm.  For  all  retrieved  positions,  the  maximum  error  in  the  lateral  positions  is  3.0  mm,  and  the  maximum  error  in  the  axial 
positions  is  1.5  mm.  The  errors  in  the  lateral  positions  increase  as  the  targets  get  closer. 
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Table  2.  Positions  of  two  targets  separated  with  different  distances 


Known  Separation 
[Ax  (mm)] 

Obj.# 

Known  Position 
[x,  v,  z  (min)] 

Retrieved  Position 
[x,  v,  z(mm)] 

Error  (mm) 

Retrieved  Separation 
[Ax  (mm)] 

11.8 

1 

27.6,  26.0,  30 

30.3,  24.4,  28.5 

2.7,  1.6,  1.5 

6.9 

2 

40.2,  26.0,  30 

37.2,  25.2,  29.5 

3.0,  0.8,  0.5 

16.8 

1 

25.1,26.0,  30 

26.4,  24.4,  28.5 

1.3,  1.6,  1.5 

14.6 

2 

42.7,  26.0,  30 

41.0,  25.2,  29.5 

1.7,  0.8,  0.5 

26.8 

1 

20.1,26.0,  30 

19.5,25.2,  29.5 

0.6,  0.8,  0.5 

27.6 

2 

47.7,  26.0,  30 

47.1,25.2,  30.5 

0.6,  0.8,  0.5 

Similar  experiments  were  also  done  when  the  target  had  a  higher  concentration  of  ink  corresponding  to  a  greater 
absorption  coefficient.  In  that  case,  smaller  errors  occurred  in  the  retrieved  position  which  should  be  due  to  higher  signal-to- 
noise  ratio. 


5.  SUMMARY  AND  DISCUSSIONS 

The  work  presented  above  demonstrate  that  TROT  is  effective  in  obtaining  three-dimensional  location  information  of  a 
single  absorptive  target  at  different  axial  locations,  and  two  absorptive  targets  separated  by  different  distances  along  a 
lateral  direction.  The  known  lateral  positions  of  the  targets  were  retrieved  within  ~2  mm  for  breast  tissue  simulating 
samples.  The  accuracy  for  assessing  the  axial  position  (depth)  depended  on  the  target  location.  High  accuracy  was 
obtained  for  target  in  the  mid-plane,  while  the  accuracy  deteriorated  as  the  target  was  moved  towards  the  source  plane  or 
the  detector  plane.  Breast  tumors  commonly  grow  further  from  the  breast  surface,  and  the  TROT  approach  is  expected 
to  be  effective  in  locating  those  even  at  an  early  stage. 

The  experiments  presented  in  this  article  were  carried  out  for  samples  in  slab  geometry.  The  approach  can  be 
adapted  to  cylindrical  geometry  as  well.  While  the  results  presented  here  are  for  absorptive  target(s),  our  preliminary 
work  show  that  TROT  will  be  effective  for  scattering  targets  as  well,  and  we  are  in  the  process  of  developing  it  for 
scattering  targets,  as  well  as,  for  targets  that  are  both  scattering  and  absorptive  in  nature.  Our  experiments  used 
continuous- wave  light,  but  the  approach  would  work  with  frequency  domain  and  time-resolved  data  as  well. 

One  important  difference  between  the  sample  used  in  our  experiments  and  human  breast  is  that  in  our  sample  the 
intervening  medium  can  be  treated  as  uniform,  while  human  breast  is  not  uniform.  The  pseudo  spectrum  is  calculated 
using  the  known  Green’s  function  of  the  background  medium.  When  the  medium  is  non-uniform,  Green’s  function  for  a 
uniform  medium  is  used  as  an  approximation.13 

The  approach  is  fast  as  no  iteration  is  involved.  The  code  used  to  perform  the  computation  is  written  in  Matlab.  With 
pre-processed  data,  a  typical  time  to  perform  the  3D  image  reconstruction  (5  cm  x  5  cm  x  6  cm)  with  a  sub-millimeter 
resolution  using  a  Pentium  2GHz  CPU  processor  and  2GB  memory,  is  approximately  10  minutes.  The  computation  time 
can  be  further  reduced  by  improving  our  algorithm,  and  using  C/C++  code. 

In  summary,  the  TROT  approach  based  on  the  concept  of  time  reversal  and  using  MUSIC  shows  promise  to  locate 
objects  in  turbid  media,  such  as  a  tumor  in  human  breast  with  useful  accuracy. 
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ABSTRACT 

Optical  imaging  using  independent  component  analysis  (OPTICA)  and  time  reversal  optical  tomography  (TROT) 
approaches  are  used  to  detect,  locate,  and  obtain  cross-section  images  of  two  tumor  pieces  inside  a  model  human  breast 
assembled  using  ex  vivo  human  breast  tissues  and  configured  as  a  semi-cylindrical  slab  of  uniform  thickness.  The 
experimental  arrangement  realized  a  multi-source  probing  scheme  to  illuminate  an  end  face  (source  plane)  of  the  slab 
sample  using  750  nm,  800  nm  and  830  nm  beams  of  laser  light.  A  multi-detector  signal  acquisition  scheme  measured 
transmitted  light  intensity  distribution  on  the  other  end  face  (detection  plane).  This  combined  multi-source  probing  and 
multi-detector  sensing  approach  culminated  in  multiple  spatial  and  angular  views  of  the  sample  necessary  for  target 
localization.  The  perturbations  in  light  intensity  distribution  in  the  detection  plane  were  analyzed  using  both  the  OPTICA 
and  the  TROT  approaches  to  obtain  locations  of  the  tumor  pieces.  A  back-projection  technique  with  OPTICA  provided 
cross-section  images  and  estimates  of  cross  section  of  the  targets  within  the  sample.  The  estimated  locations  and 
dimensions  of  targets  are  in  good  agreement  with  the  results  of  a  corroborating  magnetic  resonance  imaging  experiment 
and  known  values. 

Keywords:  Medical  and  biological  imaging,  Image  reconstruction  techniques,  Light  propagation  in  tissues,  Inverse 
problems,  Three-dimensional  image  processing,  Tomography,  Breast  cancer,  Magnetic  resonance  imaging,  Independent 
component  analysis. 


1.  INTRODUCTION 

Optical  detection  of  targets  in  a  turbid  medium  makes  use  of  the  difference  in  optical  properties,  such  as,  scattering 
coefficient,  absorption  coefficient,  index  of  refraction,  and  fluorescence  between  the  targets  of  interest  and  the 
intervening  medium  [1].  Multiple  scattering  of  light  by  the  turbid  medium  produces  a  noise  background  that  blurs  the 
image,  and  in  severe  cases  makes  direct  imaging  impossible.  Inverse  image  reconstruction  approaches  that  are 
commonly  used  to  retrieve  image  information  have  to  deal  with  the  fact  that  inverse  problems  are  ill  posed,  and  attain 
different  measures  of  success  [2]. 

We  are  developing  the  optical  tomography  using  independent  component  analysis  (OPTICA)  [3-5],  and  time  reversal 
optical  tomography  (TROT)  [6-8]  approaches  for  detecting  and  obtaining  three-dimensional  location  information  of 
target(s)  in  highly  scattering  turbid  media  in  general,  and  of  tumor(s)  in  human  breast,  in  particular.  In  this  paper  we 
present  the  results  of  our  study  of  a  “realistic  model  cancerous  breast”  formed  using  ex  vivo  human  breast  tissues  with 
two  pieces  of  tumors  embedded  within.  We  use  a  multi-source,  multi-wavelength  probing  and  multi- detector  signal 
acquisition  scheme,  and  analyze  the  resulting  data  using  both  OPTICA  and  TROT  approaches  to  obtain  images  and 
locations  of  the  tumor  pieces.  These  results  are  compared  with  magnetic  resonance  imaging  measurements  as  reference. 

2.  EXPERIMENTAL  METHODS  AND  MATERIALS 

The  experimental  arrangement  for  detecting  and  locating  tumors  in  the  realistic  model  breast  using  OPTICA  and  TROT 
is  shown  schematically  in  Fig.  1.  The  model  breast  was  assembled  using  two  pieces  of  normal  ex  vivo  female  human 
breast  tissues.  Two  pieces  of  cancerous  tissues  were  sandwiched  at  different  locations  within  the  two  normal  pieces.  The 
normal  breast  tissue  specimens  weighed  119  grams  and  127  grams  and  consisted  primarily  of  adipose  tissue,  while  each 
tumor  (infiltrating  ductal  carcinoma)  piece  weighed  approximately  1  gram.  The  sample  was  placed  inside  a  cylindrical 
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transparent  plastic  container  of  diameter  110  mm  with  a  movable  end,  which  was  moved  to  slightly  compress  the  tissue 
along  the  z-axis  and  hold  it  in  place.  The  tumor  piece  located  at  the  left  side  of  the  sample  (‘left  tumor’)  was  of  the  size 
of  8  mm  x  5  mm  x  3  mm  and  that  located  on  the  right  side  (‘right  tumor’)  was  of  the  size  10  mm  x  10  mm  x  5  mm.  The 
axial  orientation  of  the  plastic  container  and  sample  within  it  was  preserved  for  magnetic  resonance  imaging  (MRI) 
experiments  following  the  optical  measurements. 

The  optical  imaging  experiments  were  carried  out  using  a  beam  of  light  collimated  to  a  1-mm  spot  size.  Light  of  three 
different  wavelengths  750  nm,  800  nm,  and  830  nm  from  a  Ti:sapphire  laser  were  used  in  the  measurement.  The  average 
beam  power  was  maintained  at  1 0  mW  for  every  wavelength.  Multiple  source  illumination  was  realized  in  practice  by 
scanning  the  sample  in  a  26  x  16  array  of  x-y  grid  points  with  a  step  size  of  2.5  mm,  while  the  CCD  camera  imaged  the 
entire  detection  plane.  Each  illuminated  pixel  of  the  1024  x  1024  pixels  of  the  CCD  camera  could  be  regarded  as  a 
detector.  For  illumination  of  every  scanned  point  on  the  source  plane,  the  CCD  camera  recorded  the  diffusely  transmitted 
intensity  pattern  on  the  detection  plane. 

For  MRI  experiments,  the  breast  model  sample  in  the  plastic  container  was  taken  to  Memorial  Sloan-Kettering  Cancer 
Center  (MSKCC)  small  animal  MRI  facility.  The  facility  currently  utilizes  a  4.7-T  33-cm  (Bruker  BioSpin).  MR  images 
of  the  sample  were  recorded  in  2.0-mm  slice  thick  sagittal  slices. 


BC 


Ti-sapphire 
laser  oscillator 


Diode  Laser 


Figure  1.  Schematic  diagram  of  the  experimental  arrangement  (M  =  mirror,  BC=  beam  collimator,  CCD  =  charge-coupled 
device,  PC  =  personal  computer) 


3.  THEORETICAL  FORMALISM 

The  perturbation  of  the  detected  light  intensities  on  the  boundaries  of  the  medium,  the  scattered  wave  field,  due  to 
scattering  inhomogeneities  is  given,  in  Diffusion  Approximation  (DA),  by  [3-5] 

_A/ ( rd  ’  rs  ’  k)  =  j  d3r8^a  (r,  X)cG( ,  r;  X)G(r,  r  ;  l) 

+j  dr ■8D(r,  X)cV  G(?d  ,  r;  X)  ■  V  rG{  r,  r  ;  X) 

in  the  first-order  Born  approximation  assuming  that  light  diffuses  inside  the  medium.  In  Eq.(l),  Ys  and  Yd  are  the 
positions  of  the  source  and  the  detector  on  the  boundary,  respectively;  8jua  (r,  X)  =  jua  (r,  A)  -  jjao  (A)  and  8D{ r,  A)  =  D( r, 
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A)  -  D0(A)  are  the  differences  in  absorption  coefficient  and  diffusion  coefficient,  respectively,  between  the  target  at  r  and 
the  background;  c  is  the  speed  of  light  in  the  medium;  and  G( r,  r’;  A)  is  the  Green’s  function  describing  light 
propagation  from  r’  to  r  inside  the  background  turbid  medium  of  absorption  and  diffusion  coefficients  jua0  (A)  and 
D0(A )  where  A  is  the  wavelength  of  the  probing  beam.  OPTICA  Formalism  has  been  detailed  elsewhere  [3-5].  For 
TROT  algorithm,  the  experimental  data  obtained  using  the  multiple-source  and  multiple-detector  arrangement,  were  as  a 
response  matrix  K,  expressed  as [6-7] 


M  i  M  7 

K  =  {Kt  ■}=  I  Ga(Rl,Xm)rmGs(Xm,R  •)=  Zg^mSm  (2) 

m= 1  m= 1 

K  =  {Kn}=  f  Gs(Xm,Rj)rmGd(Rl,Xm)  =  f  gSmrmgdm 
J  ’  m- 1  m- 1 


where  j  =  1,2,  . . . ,  A  and  /  =  1 ,  2,  . . . .,  A  are  indices  of  sources  and  detectors,  respectively;  m=  1,2,  . . . .,  M  is  the  index 
of  the  targets  with  M  <  A;  Rp  Ri  and  Xm  are  the  positions  of  a  source,  a  detector  and  a  target  respectively;  xw  is  the 
difference  in  the  optical  parameters  (absorption  and  scattering)  of  the  target  from  that  of  the  background  medium  (the 
background  medium  may  be  uniform  or  non-uniform);  Gs(Xm,  Rj)  and  Gd(R/,  Xm)  are  the  Green’s  functions  in  the 

background  medium.  The  vector  g ^  =  [G^(XW,R1),G^(XW,R2),-  •  -,G^  (Xm,RN)]  in  Eq.  (2)  has  A  components 

[6-7].  From  K  a  time  reversal  matrix  T=  KTK  is  constructed,  and  its  eigenvalues  and  eigenvectors  are  determined. 
Leading  eigenvalues  lead  to  the  targets. 


4.  RESULTS 


4.1  OPTICA 

OPTIC A-generated  independent  intensity  distributions  on  the  detector  plane  for  the  two  tumors  using  750-nm 
probing  are  shown  in  Figure  2(a).  Similar  intensity  distributions  were  obtained  for  probing  using  light  of  other  two 
wavelengths  as  well.  The  x-y-z  locations  of  the  left  tumor  and  the  right  tumor  determined  from  fitting  these  ICA 
profiles  to  the  Green’s  function  for  all  three  wavelengths  are  shown  in  Table  I.  The  average  of  the  probed  tumor 
positions  from  all  wavelengths  is  shown  in  bold.  The  cross-section  images  of  the  two  tumors  constructed  using 
Back-projection  Fourier  transform  [5]  are  shown  in  Figure  2(b). 

Table  I.  The  coordinates  (x,  y,  z)  of  the  left  tumor  and  right  tumor 


Wavelength  (nm) 

y  (mm) 

Left  Tumor 
x  (mm) 

z  (mm) 

y  (mm) 

Right  Tumor 
x  (mm) 

z  (mm) 

750 

50.3 

54.9 

20 

26.1 

48.7 

20.2 

800 

50.5 

55.4 

17.9 

27.1 

49.4 

19.5 

830 

49.6 

54.2 

17.2 

27.1 

47.2 

21.5 

Average 

50.1 

54.8 

18.4 

26.8 

48.4 

20.4 

4.2  TROT 

The  TROT  algorithm- generated  cross-section  image  of  the  tumors  for  2,  =  750  nm  is  shown  in  Figure  3(a).  Similar 
images  were  obtained  using  other  two  wavelengths.  The  locations  were  found  to  be  (1 1.5mm,  33.1mm,  22.5  mm)  for  the 
left  tumor  and  (42.6  mm,  38.6  mm,  20.5  mm)  for  the  right  tumor.  Since  the  tumors  depth  location  is  at  the  mid  plane  - 
20-mm,  we  compare  the  horizontal  separation  between  the  two  tumors  and  find  it  to  be  35.8-mm  using  OPTICA  and 
3 1.1 -mm  using  TROT.  The  optical  results  are  in  good  agreement  with  the  MR  images  at  depth  22.5  mm  shown  3(b)  and 
3(c)  and  separation  of  35  mm. 
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x  (mm)  x  (mm) 


(b) 

Figure  2.  (a)  OPTICA  generated  intensity  distributions  pertaining  to  the  left  tumor  (LT)  and  right  tumor  (RT)  for  750-nm 
probing,  (b)  Cross-section  images  formed  by  Back-projection  Fourier  transform. 


Figure  3.  (a)  TROT  generated  cross-section  image,  (b)  Left  tumor  LT  at  slice  depth  z  =  18.5-mm,  and  (c)  Right  tumor  RT  at 
slice  depth  z  =  22. 5 -mm 
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5.  SUMMARY 


The  two  approaches  provided  comparable  information,  and  successfully  detected  and  located  the  two  tumors.  The  use  of 
multi-wavelength  method  eliminated  some  artifacts.  The  implementation  time  was  under  1  minute.  Although  we  have 
used  the  DA  in  this  case,  the  approach  can  use  other  light  propagation  models.  It  can  be  used  with  continuous  wave, 
frequency-domain  and  time-resolved  experimental  data.  Future  work  on  this  approach  will  be  directed  towards 
estimating  the  optical  properties,  as  well  as,  size  and  shape  of  targets. 
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Abstract 

Diffuse  optical  imaging  (DOI)  for  detecting  and  locating  targets  in  a  highly  scattering  turbid 
medium  is  treated  as  a  Blind  Source  Separation  (BSS)  problem.  Three  matrix  decomposition 
methods,  Independent  Component  Analysis  (ICA),  Principal  Component  Analysis  (PCA)  and 
Non-negative  Matrix  Factorization  (NMF)  were  used  to  study  the  DOI  problem.  The  efficacy  of 
resulting  approaches  was  evaluated  and  compared  using  simulated  and  experimental  data. 
Samples  used  in  the  experiments  included  Intralipid- 10%  or  Intralipid-20%  suspension  in  water 
as  the  medium  with  absorptive  or  scattering  targets  embedded. 


1.  Introduction 

Diffuse  optical  imaging  (DOI)  for  detection  and  retrieval  of  location  information  of  targets  in  a 
highly  scattering  turbid  medium  may  be  treated  as  a  blind  source  separation  (BSS)  problem  [1,  2], 
Various  matrix  decomposition  methods,  such  as,  Independent  Component  Analysis  (ICA)  [3], 
Principal  Component  Analysis  (PCA)  [4]  and  Non-negative  Matrix  Factorization  (NMF)  [5,  6] 
have  been  developed  for  solving  the  BSS  problem  and  retrieving  desired  information. 

Min  Xu  et  al.  adapted  ICA  of  information  theory  to  develop  Optical  Tomography  using 
Independent  Component  Analysis  (OPTICA)  and  demonstrated  its  application  for  diffuse 
imaging  of  absorptive,  scattering  and  fluorescent  targets  [7-11].  ICA  assumes  the  signals  from 
different  targets  to  be  independent  of  each  other,  and  optimizes  a  relevant  measure  of 
independence  to  obtain  the  ICs  associated  with  different  targets.  The  position  co-ordinates  of 
targets  in  three  dimensions  are  determined  from  the  individual  components  separately. 

PCA  assumes  that  the  PCs  contributing  to  the  signal  are  uncorrelated  and  explain  the  most 
variance  in  the  signal.  PCA  has  been  widely  used  in  various  applications,  such  as  spectroscopy 
[12],  face  recognition  [13]  and  neuroimaging  [14],  NMF  seeks  to  factorize  a  matrix  into  two 
non-negative  matrices  (component  signals  and  weights)  and  requires  the  contributions  to  signal 
and  the  weights  of  the  components  to  be  non-negative.  It  does  not  imply  any  relationship 
between  the  components.  NMF  has  also  been  widely  used  in  biological  analysis  [15],  and 
spectral  analysis  [16]. 

The  objective  of  this  study  is  to  test  and  compare  the  efficacy  of  these  algorithms  when  used 
to  solve  the  DOI  problem.  Results  are  presented  and  compared  using  simulative  data  and 
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experimental  data  using  absorptive  and  scattering  targets  embedded  in  model  scattering  media. 
Our  interest  in  solving  the  DOI  problem  derives  from  the  need  for  a  non-invasive  modality  for 
detecting,  locating,  and  diagnosing  breast  tumors  in  early  stages  of  growth. 

The  remainder  of  the  article  is  organized  as  follows.  In  Section  2,  the  formalisms  of  the  three 
methods  are  introduced.  Section  3  evaluates  the  resulting  imaging  approaches  using  simulated 
data.  The  approaches  are  further  examined  in  Section  4  for  experimental  data  acquired  using 
absorptive  and  scattering  targets  embedded  in  model  scattering  media.  Section  5  summarizes  and 
discusses  the  results. 


2.  Formalism 

2.1  Blind  Source  Separation  problem 

Blind  source  separation  (BSS),  also  known  as  blind  signal  separation  is  a  general  problem  in 
information  theory  that  seeks  to  separate  different  individual  signals  from  the  measured  signals, 
which  are  weighted  mixtures  of  those  individual  signals.  Assuming  M  individual  signals, 
Sj(t );  j  =  1,  •••  M,  are  linearly  mixed  instantaneously,  the  BSS  problem  is  modeled  as  following. 
The  dimension  of  s;  (t)  is  Ns,  the  number  of  sampling  times.  In  this  presented  study,  t  will  be 
replaced  by  spatial  positions  of  the  excitation  light  sources.  A  total  of  Nd  detectors  sense  N(/ 
different  mixtures  of  s;(t).  The  mixture  measured  by  the  ith  detector  can  be  presented  as  (t)  = 
'ZljI=1CLijSj(t),  or  X  =  AS,  in  a  matrix  notation,  where  A  6  RNdXM  is  a  mixing  or  weighting 
matrix,  5  e  RMxN\  X  E  RNdXNs ?  and  m  <  mjn  (Ns,Nd).  The  objective  of  BSS  is  to  retrieve  the 
signals  Sj  (t)  and  their  weights,  a ij.  ICA,  PC  A  and  NMF  are  statistical  analysis  methods,  used  to 
solve  the  BSS  problem. 

2.2  Diffuse  optical  imaging  problem 

In  DOI  one  measures  the  signal  at  the  sample  boundary,  which  includes  a  weighted  mixture  of 
contributions  from  embedded  targets.  One  uses  the  diffusion  approximation  [17-19]  of  the 
radiative  transfer  equation  [20,  21]  as  the  forward  model  to  describe  light  propagation  in  a  highly 
scattering  turbid  medium.  The  perturbation  in  the  light  intensity  distribution  measured  on  the 
boundary  of  the  sample  due  to  the  presence  of  the  targets  (which  are  localized  inhomogeneities 
in  the  optical  properties  within  the  sample  volume)  may  be  written,  in  the  first  order  Bom 
approximation,  as  [22,  23] 

A0(rd,rs)  =  -  f  G(rd,r)8iia(r)cG(r,rs)d3r 

-  j  SD(r)cVrG(rd,r)  ■VrG(r,rs)d3r,  (1) 

where  rs,rd,  and  r  are  the  positions  of  a  source  of  unit  power,  detector  and  target,  respectively; 
G(r,  rs)  and  G(rd,  r)  are  the  Green’s  functions  that  describe  light  propagation  from  the  source 
to  the  target  and  from  the  target  to  the  detector,  respectively;  d/na  and  SD  are  the  differences  in 
absorption  coefficient  and  diffusion  coefficient  between  the  targets  and  the  background  medium, 
respectively;  and  c  is  the  light  speed  in  the  medium. 

A  multi-source  illumination  and  multi-detector  signal  acquisition  scheme  is  used  to  acquire 
light  transmitted  through  a  scattering  medium.  For  small  absorptive  targets,  a  perturbation  data 
matrix  is  constructed  using  —A 0  for  all  sources.  The  elements  of  the  data  matrix  pertaining  to 
absorptive  targets  represented  by  the  first  term  in  Eq.  (1)  may  be  written  in  a  discrete  form  as: 
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Xtj  =  I.m=i  Gd(ri,rm)TmGs(rmtrj),  ( i  =  1,2  ,-Nd;j  =  1,2,  -,NS),  (2) 

where  r£,  r,-  and  rm  are  the  locations  of  the  r  detector,/11  source  and  /«th  target,  respectively;  A^, 
Nd  and  Mare  the  numbers  of  sources,  detectors  and  targets,  respectively;  rm  =  8na(rm)c8Vm  is 
the  optical  absorption  strength  of  the  mth  target,  8Vm  is  the  volume  of  the  target;  Gs(rm,rj )  and 
Cd(r£)  rm)  are  the  Green’s  functions  that  describe  light  propagation  from  jth  source  to  mth  target 
and  from  mth  target  to  zth  detector,  respectively.  The  number  of  targets  is  assumed  to  be  less  than 
that  of  sources  and  detectors,  M  <  min  ( Nd ,  Ns). 

The  mth  target  may  be  considered  to  be  a  virtual  source  of  strength  zmGs(rm,  ry)  excited  by 
the  real  light  source  located  at  ry .  The  data  matrix  X  =  {Xy},  may  be  considered  to  be  a  set  of 
combinations  of  light  signals  from  all  virtual  sources  mixed  by  a  mixing  matrix  { Gd(ri,rm )}. 
Therefore,  this  problem  can  be  treated  as  a  BSS  problem. 


As  the  second  term  in  Eq.  (1)  suggests,  each  scattering  target  is  represented  by  three  co- 

d 

located  virtual  sources  of  strength:  rm  dpGs(rm,rj),  where  dp  =  —  ,(p  =  x,y,z),  and  rm  = 

5D(rm)cdl^n,  is  the  optical  scattering  strength  of  the  mth  target  [8],  The  mixing  matrices  become 
[dpGd (r£,rm)},  (p  =  x,y,z ),  for  the  three  virtual  sources  generated  by  the  mth  target.  The 
elements  of  the  data  matrix  for  scattering  targets  may  be  written  as 


Xij  —  ^im=l^ip={x,y,z}  8pG  (Xi>^m)^m8pG 


(3) 


Since  one  absorptive  target  is  represented  by  one  centrosymmetric  virtual  source,  while  three 
virtual  sources  (one  centrosymmetric  and  two  dumb-bell  shaped)  represent  one  scattering  target 
[7,  8],  the  number  and  patterns  of  virtual  sources  may  be  used,  in  favorable  situations,  to  identify 
the  target  as  absorptive  or  scattering  in  nature.  In  this  paper,  only  small  targets  are  considered 
since  all  three  algorithms  are  suited  for  small  targets,  and  early  detection,  when  the  tumors  are 
more  amenable  to  treatment,  is  of  practical  interest. 


2.3  DOI  as  a  BSS  problem 

The  data  matrix  for  the  DOI  problem  may  be  written  as 


y  —  as  —  yM  A-  s  ■ 

—  no  —  Zj7M=1 

(4) 

where  A  E  RNaxM ?  5  e  rmxns^  an(j  %  ^  ^ NdxNs .  por  absorptive  targets, 

Aim  —  finrfi  Smj  —  OCmG  (v m>  ^/)? 

(5a) 

while  for  scattering  targets, 

Aim  —  Pm^pG  {Xit^m)  and  Smj  —  dmdpG  (v  m>  ^/)* 

(5b) 

{Smj}  (j  =  1,2,  ■■■ ,  Ns  )  and  {Aim}  ( i  =  1,2,  ■■■  Nd)  are  two-dimensional  intensity  distributions 
on  the  source  and  detector  planes,  respectively.  Source  and  detector  planes  are  the  boundaries  of 
the  sample  through  which  light  enters  and  exits  the  sample  volume,  respectively.  The  scaling 
factors  (3m  and  am  are  related  to  the  target  optical  strength,  rm  =  am(3m.  The  location  of  the 
target  and  the  scaling  factors  can  be  retrieved  using  a  least  squares  fitting  via 

argmin am,pm,rn&][«m  lsmj  ~  Gs(rm,rj)]2 

+  £i[An~%m  -  Gd(ri,rm)]2},  or 

(6a) 
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argminam ,/?m,rmEpGy \Pm  Smj  dpG  (j'm>  c)l 

+  -  5pGd(ri,rm)]2}},  (6b) 

for  absorptive  and  scattering  targets,  respectively.  However,  when  a  scattering  target  is  embedded 
deep  in  a  turbid  medium,  only  the  zm  dzGs(rm,rj)  virtual  source  remains  significant.  So  only 
p  =  z  may  be  used  for  fitting  in  Eq.  (6b)  [8], 

2.3.1  1C A 

OPTICA  assume  that  the  virtual  sources  are  independent  of  each  other  [8],  So  they  can  be 
retrieved  through  an  iterative  process  which  seeks  to  maximize  the  independence  among  the 
components.  In  practice,  the  independent  components  are  found  by  maximizing  some  measure  of 
non-Gaussianity,  such  as  kurtosis  (the  fourth-order  cumulant),  of  the  unmixed  components.  A 
Matlab  program  for  ICA  was  adopted  from  http://sccn.ucsd.edu/eeglab/.  The  location  of  the 
target  can  be  retrieved  by  fitting  the  independent  component  intensity  distributions  (ICIDs)  to 
Green’s  functions  or  derivatives  of  Green’s  functions  using  Eqs.  6(a)  and  6(b). 

2.3.2  PC  A 


PCA  assumes  that  the  virtual  sources  are  uncorrelated  so  that  the  correlation  (covariance) 
between  them  is  ideally  zero,  and  minimal  in  practice.  The  covariance  matrix  of  S,  cov(S ) 
should  be  diagonal.  The  general  process  of  PCA  is  as  follows.  The  data  matrix  X  =  AS  +  N, 
where  N  is  random  noise  added  to  the  data,  A  and  S  the  same  as  defined  in  Eq.  (4).  When  S  is 
mean  centered,  elements  of  the  mean-centered  matrix  S'  are  defined  as 


cf  _  c  _  ^  ^ 

°  mj  ~  J mj  N 


mj- 


Similarly 


Y'  —  Y  1  vN>  V 

A  ij  —  Aij  N<j= 1  A‘/' 


(7a) 

(7e) 


PCA  looks  for  a  matrix  P  that  decomposes  X  into  virtual  sources,  S  =  PX.  It  also  holds  that 
S'  =  PX',  since  P  is  just  a  rotation  matrix  which  does  not  change  the  center  of  the  data. 

cov(S )  =  S’S,T  =  C PX')(PX')T  =  PX'X,TPT  =  A,  (8) 


where  A  =  diag{A1,A2,-‘4}  ■  The  eigenvalues  Xm  are  variances  in  the  covariance  matrix. 
Therefore,  X'X'T PT  =  PT A  ,  where  PT  is  orthonormal.  PCA  is  realized  by  eigenvalue 
decomposition  (EVD)  of  the  covariance  matrix  of  X.  The  eigenvectors  with  leading  eigenvalues 
(largest  variances)  are  selected  to  be  the  PCs  using  the  L-curve  [24], 

Since,  X  =  PTS  ~  AS,  A  is  determined  as  a  matrix  including  only  PCs.  S  is  calculated  as 
S  «  (AT  A)~1AtX  .  Rows  of  S  and  columns  of  A  represent  principal  component  intensity 
distributions  (PCIDs)  on  the  source  plane  and  detector  plane,  respectively,  and  are  proportional 
to  the  images  of  the  virtual  sources  projected  on  the  source  and  detector  planes.  The  target 
positions  are  determined  using  Eq.  (6). 


2.3.3  NMF 


NMF  is  a  group  of  multivariate  analysis  algorithms  that  factorize  a  matrix  X  into  A  and  S:  X=  AS, 
A  and  S  are  non-negative  [6].  Unlike  ICA  and  PCA,  NMF  does  not  imply  any  relationship 
between  the  retrieved  components;  instead,  it  just  enforces  non-negativity  constraints  on  A  and  S. 
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There  are  various  algorithms  developed  to  solve  NMF,  such  as  the  multiplicative  update  method 
[5]  and  alternating  least  squares  method  [25,  26], 

implementation  of  NMF,  A  and  S  can  be  found  by  minimizing 
\\X  —  AS ||2  as  the  cost  function,  where  A  >  0  and  S  >  0,  using 

(9a) 
(9b) 

The  alternating  least  squares  implementation  of  NMF  uses  alternate  least  squares  steps  to 
estimate  A  (or  S ),  and  use  that  estimate  to  optimize  S  (or  A),  and  keep  repeating  the  alternative 
steps  until  the  desired  optimization  is  obtained.  Non-negativity  is  ensured  by  setting  any  negative 
element  of  A  or  S  equal  to  0. 

A  NMF  toolbox  was  obtained  from  http://cogsvs.imm.dtu.dk/toolbox/  to  perform  NMF 
computation.  A  built-in  command  nnmf  is  also  available  in  Matlab  (R201  la). 

NMF  algorithm  requires  that  the  non-negativity  assumption  must  hold  in  the  problem.  In 
particular,  for  absorptive  targets,  when  X  is  constructed  with  —A tp,  rm  should  be  positive,  i.e.  the 
targets  should  be  more  absorbing  than  the  background.  If  the  targets  have  weaker  attenuation 
properties  than  the  background,  X  needs  to  be  constructed  with  +A0  instead.  For  scattering 
targets,  X  should  be  treated  similarly  to  keep  its  elements  positive. 

When  NMF  is  applied  to  a  scattering  target,  only  the  centrosymmetric  component  shows  up 
properly,  since  the  other  two  components  have  dumb-bell  shape  which  includes  negative  values 
[8],  So  without  any  prior  knowledge  or  some  other  experimental  means  to  assess  if  the  target  is 
absorptive  or  scattering,  NMF  may  not  distinguish  between  the  two  possibilities. 

The  decomposition  methods  can  be  applied  with  different  sample  geometries  such  as  slab  and 
cydrindrical  geometries,  and  different  measurement  domains  such  as  time -resolved  domain, 
frequency  domain  and  continuous  wave  (CW).  In  this  article,  Green’s  functions  for  slab 
geometry  [23]  with  CW  measurement  were  used  for  simulation  and  experiments. 


In  the  multiplicative  update 
the  square  of  Euclidean  distance 
the  multiplicative  update  rule 


Aik  <-  A; 


ik 


Q(£}ik 

C asst ) 


ik 


kj  kj  (ATAS)k]- 


3.  Simulation 

The  sample  was  considered  to  be  a  40-mm  thick  uniform  scattering  slab  with  lateral  dimension 
of  80  mm  x  80  mm,  as  shown  in  Fig.  1.  Its  absorption  and  diffusion  coefficients  were  taken  to  be 

fia  =  0.003  mm'1  and  D  =  -  mm  (transport  mean  free  path,  lt=  1  mm),  respectively,  which  are 

similar  to  the  average  value  of  those  parameters  for  human  breast  tissue.  An  absorptive  and  a 
scattering  point  target  were  placed  at  (50,  60,  15)  mm  and  (30,  30,  25)  mm,  respectively.  The 
index  of  refraction  n  of  the  medium  was  taken  to  be  1.33.  The  speed  of  light  is  2.998  x  108  m/s, 
or  299.8  mm/ns  in  vacuum,  and  225.4  mm/ns  in  the  medium.  The  absorption  coefficient  of  the 
absorptive  target  was  set  to  be  higher  than  the  background  with  A/ra  =  0.001  mm'1,  while  the 
diffusion  coefficient  was  taken  to  be  the  same  as  that  of  background.  The  diffusion  coefficient  of 
the  scattering  target  was  set  to  be  lower  than  the  background  (higher  scattering  coefficient)  with 
A D  =  -0.1  mm  (/,  =  0.7  mm),  while  the  absorption  coefficient  was  taken  to  be  the  same  as  the 
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background.  The  volumes  of  both  targets  are  set  to  be  8  mm  .  The  optical  strengths  of  the 
absorptive  and  scattering  targets  were  A/iacAV  =  1.803  mm  /ns,  and  ADcAV  =  -180.3  mm  /ns, 
respectively.  The  incident  CW  beam  step  scanned  the  sample  at  21  x  21  grid  points  covering  an 
80  x  80  mm2  area,  with  a  step  size  of  4  mm.  Light  on  the  opposite  side  was  recorded  at  41  x  41 
grid  points  covering  the  same  area.  Multiplicative  Gaussian  noise  of  5%  was  added  to  the 
simulated  data.  The  data  matrix  X  was  then  obtained  using  Eq.  (2)  directly,  and  analyzed  using 
the  three  different  algorithms. 


Fig.  1 :  Light  intensity  distribution  on  the  detector  plane  is  recorded  when  a  point  source  scans  on  the  source  plane. 


3.1  ICA  Analysis 

One  independent  component  for  the  absorptive  target  and  three  independent  components  for  the 
scattering  target  were  retrieved  by  ICA.  The  independent  component  intensity  distributions 
(ICIDs)  on  the  detector  plane  are  shown  in  Figs.  2(a),  2(c),  2(d),  and  2(e).  Similar  ICIDs  were 
obtained  on  the  source  plane.  Fig.  2(g)  shows  the  centrosymmetric  ICID  for  the  scattering  target, 
and  Fig.  2(i)  shows  the  ICID  for  the  absorbing  target. 

The  components  in  either  the  detector  plane  or  the  source  plane  can,  in  principle,  be  used  to 
extract  position  and  optical  strength  of  the  target(s).  However,  in  our  experimental  arrangement 
signal  is  collected  by  a  1024  x  1024  pixels  CCD  camera,  while  the  source  plane  is  scanned  in  an 
x-y  array  of  points,  which  is  much  smaller  than  the  number  of  pixels  in  the  CCD  camera. 
Consequently,  the  resolution  in  the  detector  plane  is  much  better,  and  the  data  set  more  robust 
than  the  source  side.  So,  we  used  the  images  on  the  detector  plane  for  retrieving  target 
information  using  experimental  data.  While  it  would  not  matter  in  simulation,  to  be  consistent 
with  experimental  situations,  we  employed  detector  plane  images  when  using  simulated  data  as 
well  for  all  three  algorithms.  Table  I  lists  the  locations  and  strengths  of  the  absorptive  and 
scattering  targets  retrieved  by  fitting  the  spatial  intensity  profile  of  the  centrosymmetric 
components  on  the  detector  plane  to  Green’s  functions  and  derivatives  of  Green’s  functions 
using  Eq.  6(a)  and  Eq.  6(b),  respectively,  as  shown  in  Fig.  2(b)  and  Fig.  2(f).  Fig.  2(h)  and  Fig. 
2(j)  show  the  corresponding  fits  to  the  profiles  on  the  source  plane. 
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Detector  Plane  Images  and  Profiles 


20 
|  40 
^60 
80 
(a) 


(c)  x(mm) 


20  40  60  80 


20 

f  40 

^60 

80 

_ 

20  40  60  80 
(d)  x  fmm) 


x(mrn) 

Source  Plane  Images  and  Profiles 


Fig.  2:  ICA  extracted  two-dimensional  intensity  distribution  on  the  detector  plane  of:  (a)  the  centrosymmetric 
component,  (c),  (d)  dumb-bell  shaped  components  of  the  scattering  target;  and  (e)  the  absorptive  target.  Similar 
intensity  distribution  on  the  source  plane  of:  (g)  the  centrosymmetric  component  of  the  scattering  target,  and  (i)  the 
absorptive  target  for  comparison.  Fits  to  the  spatial  intensity  profile  on  the  detector  plane  along  the  white  dashed  line 
(shown  in  figures)  of  the  centrosymmetric  component  of  the  scattering  target  is  shown  in  (b),  and  that  of  the 
absorptive  target  is  shown  in  (f).  Corresponding  fits  to  spatial  profiles  on  the  source  plane  are  displayed  in  (h)  and 
(j),  respectively. 


3.2  PCA  Analysis 

Eigenvalue  equation  of  the  covariance  matrix  of  X  was  solved.  The  eigenvalues  found  by  PCA 
were  sorted  in  descending  order.  Fig.  3  shows  a  plot  of  leading  20  eigenvalues  on  a  logarithmic 
scale. 
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Fig.  3:  A  logarithmic  plot  of  the  first  20  PC  A  eigenvalues 

First  four  leading  eigenvalues  were  selected  for  PCs.  The  corresponding  PCIDs  were 
calculated.  The  PCIDs  on  the  detector  plane  are  shown  in  Fig.  4.  Similar  images  for  PCIDs  on 
the  source  plane  were  obtained.  The  scattering  target  has  one  centrosymmetric  (Fig.  4(a)) 
component  and  two  dumb-bell  shaped  (Fig.  4(c)  and  Fig.  4(d))  components,  while  the  absorptive 
target  has  only  one  component  (Fig.  4(e)). 
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Fig.  4:  PCA-extracted  two  dimensional  intensity  distribution  on  the  detector  plane  of:  (a)  the  centrosymmetric 
component,  and  (c),  (d)  dumb-bell  shaped  components  of  the  scattering  target;  (e)  the  absorptive  target.  Green’s 
function  fits  to  the  spatial  intensity  profiles  along  the  dashed  line  (shown  in  figures)  of  the  (b)  centrosymmetric 
component  of  the  scattering  target  and  (f)  absorptive  target,  respectively,  to  retrieve  the  locations  of  the  two  targets. 

Fig.  4(b)  and  Fig.  4(f)  show  fits  to  the  spatial  intensity  profile  of  the  centrosymmetric 
component  of  the  scattering  target  and  that  of  the  absorptive  target,  respectively,  to  retrieve  the 
locations  of  the  two  targets.  The  locations  and  optical  strengths  of  the  targets  retrieved  by  PCA 
are  also  shown  in  Table  I. 

3.3  NMF  Analysis 

The  mixing  matrix  and  virtual  sources  were  retrieved  from  the  data  matrix  X  using  NMF  as 
explained  in  Section  2.3.3.  As  in  the  other  two  approaches,  only  one  component  is  extracted  for 
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the  absorptive  target.  Since  NMF  has  a  non-negativity  constraint,  only  the  centrosymmetric 
component  for  the  scattering  target  is  obtained.  Non-negative  component  intensity  distributions 
(NCIDs)  on  detector  planes  are  shown  in  Fig.  5.  Similar  images  for  NCIDs  on  source  plane  were 
also  obtained  using  the  virtual  sources  in  S.  The  results  are  also  shown  in  Table  I. 


Fig.  5:  NMF-extracted  two  dimensional  intensity  distribution  on  the  detector  plane  of:  (a)  the  centrosymmetric 
component  of  the  scattering  target;  (c)  the  absorptive  target.  Fits  to  the  corresponding  spatial  intensity  profiles  along 
the  dashed  line  (shown  in  figures)  are  shown  in  (b)  and  (d),  respectively. 

3.4  Results  and  discussion 

The  positions  and  optical  strengths  of  the  targets  retrieved  by  ICA,  PCA  and  NMF  algorithms  are 
shown  in  Table  I,  and  compared  to  the  known  values.  The  retrieved  results  using  all  three 
algorithms  from  this  simulated  data  are  in  excellent  agreement  with  the  known  values. 


Table  I.  Positions  and  optical  strengths  retrieved  using  simulated  data  and  ICA,  PCA  and  NMF 

algorithms 


Target 

Known 

position 

(mm) 

Algorithm 

Fitted 

position 

(mm) 

Error 

(mm) 

Known* 

strength 

Fitted* 

strength 

Error 

(%) 

Sea. 

(30,  30,  25) 

ICA 

(29.9,  30.0,  25.1) 

(0.1,0,  0.1) 

-180.3 

-179.9 

0.22 

PCA 

(30.0,  30.0,  25.0) 

(0, 0,  0) 

-180.3 

-180.1 

0.11 

NMF 

(30.0,  30.0,  25.0) 

(0, 0,  0) 

-180.3 

-178.5 

1 

Abs. 

(50,  60,  15) 

ICA 

(50.1,60.2,  15.0) 

(0.1,  0.2,  0) 

1.803 

1.826 

1.28 

PCA 

(50.1,60.1,  14.9) 

(0.1,  0.1,  0.1) 

1.803 

1.812 

0.5 

NMF 

(50.1,60.1,  15.0) 

(0.1,  0.1,0) 

1.803 

1.803 

0 

*  The  unit  for  absorption  strength  of  the  target  is  mm3/ns  and  for  scattering  strength  is  mm5/ns. 


9 


4.  Experiments 

4.1  Experimental  materials  and  methods 

In  this  Section,  the  algorithms  are  evaluated  using  experimental  data  for  absorptive  and  scattering 
targets  embedded  in  model  scattering  media  whose  absorption  and  scattering  properties  are 
adjusted  to  mimic  the  average  values  of  those  parameters  for  human  breast  tissues.  Two  different 
experiments  were  carried  out  with  two  different  samples.  The  first  sample  used  a  250  mm  x  250 
mm  x  50  mm  transparent  plastic  container  filled  with  Intralipid- 10%  suspension  in  water  as  the 
background  medium.  The  concentration  of  Intralipid- 10%  was  adjusted  to  provide  [27,  28]  an 
absorption  coefficient  of  jua  ~  0.003  mm'1,  and  a  transport  mean  free  path  lt  ~  1.43  mm  at  785  nm. 
The  second  sample  used  a  similar  container  with  dimension  of  250  mm  x  250  mm  x  60  mm  filled 
with  Intralipid-20%  suspension  in  water.  The  concentration  of  Intralipid-20%  was  adjusted  to 
provide  [27,  28]  jua  ~  0.003  mm'1,  and  /,  ~  1  mm  at  785  nm.  These  optical  parameters  of  the 
medium  were  selected  to  be  similar  to  the  average  values  of  those  parameters  for  human  breast 
tissue.  The  thickness  of  the  samples  was  also  comparable  to  that  of  a  typical  compressed  female 
human  breast. 

In  the  first  experiment,  two  absorptive  targets  were  embedded  in  the  medium.  The  targets  were 
~  10-mm  diameter  glass  spheres  filled  Indocyanine  green  (ICG)  dye  dissolved  in  Intralipid-20% 
suspension  in  water  to  obtain  an  absorption  coefficient  jua  =  1.15  mm'1  at  785  nm,  and  to  match  the 
background  scattering  coefficient  of  2.11  mm'1.  The  targets  were  placed  at  (57.2,  18.1,  20.0)  mm 
and  (19.9,  48.1,  25.0)  mm,  respectively. 

In  the  second  experiment,  two  scattering  targets  were  embedded,  which  were  also  ~  10  mm 
diameter  glass  spheres,  filled  with  Intralipid-20%  suspension  in  water.  The  transport  mean  free 
path,  lt  was  adjusted  to  be  0.25  mm,  with  scattering  coefficient  jus  ~  1 1  mm'1,  and  absorption 
coefficient  jua  same  as  the  background  medium.  The  targets  were  placed  in  the  middle  plane  (z  = 
30  mm)  in  the  container  with  a  lateral  distance  of  40  mm  from  each  other  (center  to  center). 

The  experimental  setup  is  schematically  shown  in  Fig.  6.  A  10-mW  785-nm  diode  laser  beam 
was  used  to  illuminate  the  first  sample,  while  a  100-mW  790-nm  diode  laser  beam  was  used  for  the 
second  sample.  The  input  surface  (source  plane)  of  the  samples  was  scanned  across  the  laser  beam 
in  an  x-y  array  of  grid  points  to  realize  the  multi-source  interrogation  of  the  samples.  The 
transmitted  light  from  the  exit  surface  (detector  plane)  was  recorded  by  a  1024  pixel  xl024  pixel 
(pixel  size  =  24  \m)  CCD  camera  (Photometries  CH350)  equipped  with  a  60-mm  focal-length 
camera  lens.  Each  pixel  of  the  CCD  camera  can  be  considered  to  be  a  detector  implementing  the 
multi-detector  signal  acquisition  arrangement.  A  set  of  16-bit  1024  pixel  xl024  pixel  images  were 
acquired.  The  two  samples  were  scanned  in  an  array  of  1 1  x  12  and  11x15  grid  points,  respectively, 
with  a  step  size  of  5  mm  in  both  cases.  The  processes  of  scanning  and  data  acquisition  were 
controlled  by  a  personal  computer.  At  all  scan  positions,  raw  transillumination  images  of  the 
samples  were  recorded  by  the  computer  for  further  analysis. 
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Fig.  6:  A  schematic  diagram  of  the  experimental  arrangement  used  for  imaging  objects  embedded  in  a  turbid 
medium.  The  inset  at  the  bottom  shows  the  2D  array  in  the  input  plane  that  was  scanned  across  the  incident  laser 
beam;  and  the  inset  to  the  right  shows  a  typical  raw  image  recorded  by  the  CCD.  (CCD  =  charge  coupled  device,  PC 
=  personal  computer) 


4.2  Analysis  and  Results 

A  region  of  interest  (ROI)  was  cropped  out  from  each  image.  Then  every  5x5  pixels  in  each 
cropped  image  were  binned  to  one  pixel  to  enhance  signal-to-noise  ratio.  A  background  image  was 
generated  by  calculating  an  average  image  for  all  scan  positions  to  approximate  the 
transillumination  image  without  target(s)  embedded. 

This  averaging  method  for  generating  background  image  is  suitable  for  small  targets  used  in 
our  experiments,  as  the  ratio  of  the  volume  of  the  sample  to  that  of  the  target  was  quite  high 
(~500:1).  For  in  vivo  imaging  of  tumors  in  early  stages  of  growth,  the  breast-to-tumor  volume 
ratio  will  be  similarly  high  and  the  averaging  method  will  be  applicable.  Alternative  approaches 
for  generating  a  background  image  include  using  image  of  (a)  a  phantom  that  has  the  same 
average  optical  properties  as  the  sample  [29];  (b)  the  healthy  contralateral  breast  for  breast 
imaging  [30];  and  (c)  the  sample  obtained  using  light  of  wavelength  for  which  the  target(s)  and 
the  background  have  identical  optical  properties  [31].  Still  another  approach  is  to  compute  the 
background  using  an  appropriate  forward  model  [32],  A  more  detailed  discussion  of  this 
important  issue  appears  in  one  of  our  earlier  publications  [33], 

The  background  image  was  also  cropped  and  binned  corresponding  to  the  ROI  for  each  scan 
position.  Perturbation  in  the  light  intensity  distribution,  A(p  due  to  targets  in  each  image  was  found 
by  subtracting  the  background  image  from  the  image.  The  data  matrix  X  was  then  constructed 
using  the  light  intensity  perturbations  at  all  scan  positions.  ICA,  PCA,  and  NMF  decomposition 
algorithms  were  performed  on  the  data  matrix  separately.  Results  are  shown  and  discussed  below. 

4.2.1  Absorptive  targets 

The  images  on  the  detector  plane  obtained  using  the  ICA,  PCA,  and  NMF  algorithms  are  shown 
in  Fig.  7,  Fig.  8,  and  Fig.  9,  respectively.  Similar  images  on  the  source  plane  were  also  obtained 
using  all  three  algorithms.  The  right  side  of  each  figure  shows  the  corresponding  spatial  intensity 
profile.  Locations  of  the  targets  are  extracted  from  fits  to  these  spatial  intensity  profiles,  as 
described  in  Section  2.3  using  Eq.  (6).  The  results  are  presented  in  Table  II.  In  Fig.  7,  images  on 
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the  source  plane  are  shown  in  (e)  and  (g),  and  Green’s  function  fits  to  their  spatial  profiles  are 
shown  in  (f)  and  (h)  for  comparison. 


Detector  Plane  Images  and  Profiles 
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Source  Plane 
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0  20  40  60 

(e)  x  (mm) 


Images  and  Profiles 


(g)  x  (mm) 


Fig.  7:  ICA-generated  ICIDs  on  the  detector  plane  are  shown  in  (a)  and  (c);  corresponding  Green’s  function  fits  to 
the  horizontal  spatial  profiles  through  the  dashed  lines  are  shown  in  (b)  and  (d).  ICIDs  on  the  source  plane  are 
shown  in  (e)  and  (g);  corresponding  Green’s  function  fits  to  the  horizontal  spatial  profiles  through  the  dashed  line 
are  shown  in  (f)  and  (h). 
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Fig.  8:  PCIDs  on  the  detector  plane  are  shown  in  (a)  and 
horizontal  spatial  profiles  through  the  dashed  line  are  shown  ii 
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(c);  and  corresponding  Green’s  function  fits  to  the 
(b)  and  (d). 
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Fig.  9:  NCIDs  on  the  detector  plane  are  shown  in  (a)  and  (c);  corresponding  Green’s  function  fits  to  the  horizontal 
spatial  profiles  through  the  dashed  line  are  shown  in  (b)  and  (d). 


It  follows  from  the  comparison  of  in  Table  II  that  the  positions  retrieved  by  all  three 
algorithms  are  in  good  agreement  with  the  known  positions.  The  errors  in  the  retrieved  locations 
(x,  y,  z)  of  the  two  targets  were  within  1.7  mm.  The  PCIDs  were  not  totally  separated.  Some 
“residue”  was  observed  in  one  PCID  from  the  other.  ICA  and  NMF  separated  two  components 
from  this  dataset  more  clearly. 
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Table  II.  Known  positions  vs.  retrieved  positions  of  the  absorptive  targets  using  ICA,  PCA  and  NMF 

algorithms 


Target 

Known  position 

(mm) 

Algorithm 

Fitted  position  (mm) 

Error  (mm) 

(57.2,  18.1,20) 

ICA 

(57.4,  18.2,21.5) 

(0.2,  0.1,  1.5) 

1 

PCA 

(57.4,  18.2,  20.6) 

(0.2,  0.1,  0.6) 

NMF 

(57.4,  18.2,  19.5) 

(0.2,  0.1,  0.5) 

(19.9,  48.1,25) 

ICA 

(18.2,  46.7,  24.7) 

(1.7,  1.4,  0.3) 

2 

PCA 

(18.2,  47.6,  25.9) 

(1.7,  0.5,  0.9) 

NMF 

(18.2,  47.6,  23.3) 

(1.7,  0.5,  1.7) 

4.2.2  Scattering  targets 


The  “images”  corresponding  to  the  centrosymmetric  components  of  the  virtual  sources  (targets) 
on  the  detector  plane  obtained  using  the  ICA,  PCA,  and  NMF  algorithms  are  shown  in  Fig.  10, 
Fig.  11,  and  Fig.  12,  respectively.  Similar  images  on  the  source  plane  were  also  obtained.  The 
right  side  of  each  figure  shows  the  corresponding  spatial  intensity  profile.  Locations  of  the 
targets  are  extracted  from  fits  to  these  spatial  intensity  profiles,  as  described  in  Section  2.3  using 
Eq.  (6).  The  results  are  presented  in  Table  III. 


Fig.  10:  ICA-generated  ICIDs  on  the  detector  plane  are  shown  in  (a)  and  (c);  corresponding  Green’s  function  fits  to 
the  horizontal  spatial  profiles  through  the  dashed  line  are  shown  in  (b)  and  (d). 
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Fig.  11:  PCIDs  on  the  detector  plane  are  shown  in  (a)  a 
spatial  profiles  through  the  dashed  line  are  shown  in  (b) 


(c);  corresponding  Green’s  function  fits  to  the  horizontal 

i  (d). 


Fig.  12:  NCIDs  on  the  detector  plane  are  shown  in  (a)  and  (c);  corresponding  Green’s  function  fits  to  the  horizontal 
spatial  profiles  through  the  dashed  line  are  shown  in  (b)  and  (d). 

Both  targets  were  detected  by  all  three  algorithms.  The  target  locations  retrieved  by  three 
algorithms  are  shown  in  Table  III,  and  compared  with  known  locations.  Overall,  all  three 
algorithms  detect  and  locate  the  scattering  and  the  absorptive  targets  with  good  accuracy,  the 
maximum  deviation  of  any  one  coordinate  from  the  known  value  being  ~3  mm.  Since  the 
maximum  difference  between  the  known  and  retrieved  position  coordinates  was  larger  for  the 
scattering  targets,  we  calculated  the  squared  correlation  coefficient  y  to  assess  the  fitting  quality. 
NMF  retrieves  the  position  coordinates  better  (within  0.5  mm)  for  the  scattering  target  #2  than 
done  by  ICA  and  PCA  (deviation  from  known  values  being  between  2-3  mm).  NMF  retrieved 
the  position  coordinates  for  target  #1  with  3.0  mm  error  in  z  direction,  which  is  not  as  good  as 
that  done  by  ICA  and  PCA.  But  y  is  0.783  and  0.778  in  the  fittings  for  ICA  and  PCA, 
respectively,  as  compared  to  0.993  for  NMF,  indicating  that  the  quality  of  the  fitting  is  better  for 
NMF.  The  quality  of  fitting  is  presumably  affected  by  the  efficacy  of  decomposition.  The 
decomposed  NCIDs  by  NMF  were  more  “clean”  than  those  decomposed  by  ICA  and  PCA.  We 
ascribe  the  observed  higher  errors  in  ICA  and  PCA  estimates  of  the  position  coordinates  of  the 
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scattering  target  #2  than  the  NMF  estimates  to  the  interference  from  the  other  virtual  source 
(corresponding  to  target  #1)  in  ICA  (Fig.  10(c))  and  PCA  (Fig.  11(c))  images.  It  is  commonly 
believed  that  errors  in  locating  a  scattering  target  are  higher  than  that  for  locating  an  absorptive 
target,  and  the  results  of  this  study  conform  to  that  notion. 


Table  III:  Known  positions  vs.  retrieved  positions  of  the  scattering  targets  using  ICA,  PCA  and  NMF 

algorithms 


Target 

# 

Known  position 

(mm) 

Algorithm 

Fitted  position 

(mm) 

Error  (mm) 

(13.0,  28.0,  30.0) 

ICA 

(12.6,  28.7,  29.1) 

(0.4,  0.7,  0.9) 

1 

PCA 

(12.6,  28.7,  28.6) 

(0.4,  0.7,  1.4) 

NMF 

(12.0,  28.5,33.0) 

(1.0,  0.5,  3.0) 

(53.3,  28.5,  30.0) 

ICA 

(51.0,31.8,  26.8) 

(2.3,  3.3,  3.2) 

2 

PCA 

(50.9,31.8,26.7) 

(2.4,  3.3,  3.3) 

NMF 

(53.3,28.0,  30.3) 

(0.0,  0.5,  0.3) 

5.  Summary  and  Discussion 

Diffusive  optical  imaging  was  modeled  as  a  BSS  problem.  ICA,  PCA  and  NMF  were  used  to 
decompose  the  data  matrix,  and  locate  the  targets  embedded  in  a  highly  scattering  turbid 
medium.  Only  the  components  corresponding  to  the  targets  were  extracted  from  a  large  dataset 
for  target  detection  and  localization. 

It  may  be  instructive  to  compare  the  objectives,  scope  and  computational  complexity  of  these 
decomposition  methods  with  model-based  reconstruction  methods.  Decomposition  methods 
obtain  the  3-D  locations  of  targets  (the  number  of  targets  are  generally  small).  Based  on  the 
retrieved  locations,  the  methods  may  then  be  further  extended  to  retrieve  size  and  optical 
property  information  of  the  targets  [9].  The  common  practice  of  model-based  inverse 
reconstruction  methods  is  to  discretize  the  sample  volume  into  N  x  N  x  N  voxels,  and  estimate 
absorption  and/or  scattering  coefficient  in  each  voxel  iteratively.  Voxels  with  significantly 
different  optical  properties  than  the  surrounding  are  regions  of  interest,  and  may  be  identified  as 
targets.  While  estimating  the  optical  properties,  the  forward  model  is  solved  repeatedly  to 
calculate  the  intensity  of  the  multiply-scattered  light  on  the  sample  boundary.  The  difference 
between  the  intensity  of  the  multiply  scattered  light  predicted  by  the  forward  model  and  the 
experimental  measurements  is  minimized  by  seeking  an  optimal  set  of  the  optical  properties  of 
every  voxel  in  the  sample  volume.  The  number  of  variables  thus  is  on  the  order  A3.  To  determine 
location(s)  of  target(s)  in  three  dimensions,  the  decomposition  methods  process  the  data  matrix 
to  retrieve  the  main  components  (A  and  S).  Here  A  and  S  are  two-dimensional  matrices  with  the 
number  of  unknowns  on  the  order  of  N  .  The  number  of  unknowns  is,  hence,  reduced  N  times  in 
the  decomposition  methods  compared  to  the  model-based  approaches,  which  leads  to  a 
substantial  saving  in  the  computational  time  when  N  is  large.  No  repeated  solution  of  the  forward 
model  is  involved  in  decomposition  methods.  Consequently,  decomposition  methods  are 
considerably  faster. 
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A  comparison  of  the  computational  complexity  of  these  two  types  of  approaches  may  shed 
further  light  on  their  relative  computation  economy.  For  a  model-based  iterative  reconstruction 
method,  an  equation  of  the  form  b  =  Wx  is  solved  to  find  the  targets,  where  W  is  a  weight  matrix 
of  size  NdNs  x  Nv,  Nj,  Ns,  and  Nv  are  the  numbers  of  detectors,  sources  and  voxels,  respectively,  b 
is  an  NdNs  x  1  vector  describing  the  perturbation  in  the  detected  light  intensity  due  to  the 
presence  of  targets,  and  x  is  the  perturbation  in  the  optical  properties  from  the  background  values 
with  dimension  of  Nv  x  1 .  The  computational  complexity  is  typically  0(NdNsNv2)  for  a  single 
iteration.  For  the  decomposition  approach,  b  is  written  as  a  2-D  matrix  X  with  dimension  Nd  x  Ns. 
To  decompose  matrix  X,  the  computational  complexity  per  iteration  is  typically  of  order  0{NdNu ) 
for  ICA  [34],  and  0(NdNsNk)  for  NMF  [35],  where  Nk  is  the  number  of  components  that  relates  to 
the  number  of  targets  and  is  usually  a  small  number.  For  PCA  using  SVD,  the  complexity  is 
0(Ns2Nk)  [34],  The  computational  complexity  of  the  intrinsic  iterative  process  involved  in  the 
matrix  decomposition  algorithms  is  much  lower  than  that  in  the  model-based  inverse 
reconstruction  methods. 

All  three  matrix  decomposition  methods  presented  in  this  manuscript  can  potentially  be  used 
in  in-vivo  real-time  breast  cancer  imaging.  The  three  algorithms  have  different  assumptions, 
which  may  lead  to  different  favored  conditions.  In  this  study,  the  algorithms  were  evaluated 
using  simulative  and  experimental  data  using  model  scattering  media  and  absorptive  and 
scattering  targets.  The  (x,  y,  z)  positions  of  the  targets  were  retrieved  with  good  accuracy.  The 
decomposition  provided  by  ICA  is  “cleaner”  than  that  of  the  PCA.  PCA  did  not  clearly  separate 
the  two  absorptive  targets  used  in  the  first  experiment.  NMF  decomposition  seems  to  provide 
residue-free  “cleaner”  images  than  the  other  two  methods  in  this  study.  However,  since  NMF  is 
based  on  non-negativity  assumption,  the  results  might  deteriorate  when  such  a  non-negativity 
assumption  does  not  hold  well.  While  continuous  wave  measurements  were  used  in  the  work 
presented  in  this  article,  the  approaches  could  be  used  with  frequency  domain  and  time-domain 
measurements  as  well. 

The  work  presented  here  focuses  on  detecting  and  locating  small  targets,  which  derive 
impetus  from  the  need  to  detect  tumors  in  early  stages  of  growth  when  those  are  more  amenable 
to  treatment.  All  three  methods  are  applicable  for  extended  targets  as  well,  and  are  expected  to 
provide  the  “center  of  optical  strength”  as  the  location  of  the  target. 

All  three  approaches  are  applicable  for  both  scattering  and  absorbing  targets,  and  may  be 
used  in  clinical  setting.  The  contrast  between  a  tumor  and  surrounding  normal  tissue  can  be  due 
to  differences  in  absorption,  scattering,  or  both  absorption  and  scattering  properties  and  may 
depend  significantly  on  the  wavelength  of  light  used.  However,  a  priori  knowledge  of  the  optical 
characteristics  (absorptive  or  scattering)  is  not  crucial.  As  has  been  shown  [Eq.  (2)  and  Eq.  (3)] 
the  expression  for  elements  of  the  data  matrix  for  absorptive  targets  involves  Green’s  Functions 

G,  while  that  for  scattering  targets  involves  dG/dz  ~  —kG^  where  k  =  ^ pa/D  in  CW  [9].  This 
relationship  with  G  provides  basis  for  detection  and  localization  of  target(s),  whether  contrast  is 
due  to  absorption,  scattering,  or  both.  We  are  using  transillumination  geometry,  which  is  one  of 
the  approaches  used  by  other  researchers,  and  adequate  signal  for  in  vivo  breast  imaging  is 
obtained  [29,  36-39], 

In  this  article,  we  presented  results  when  the  approaches  were  used  to  detect  and  obtain  three- 
dimensional  location  information  of  the  targets.  We  have  demonstrated,  while  developing 
OPTICA  [11]  that  a  back-projection  formalism  can  be  further  implemented  to  get  a  cross-section 
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image  of  the  target  [1 1],  or  the  retrieved  target  locations  can  be  fed  into  other  DOI  methods  as  a 
priori  information  to  get  three-dimensional  tomographic  images.  Since  the  approaches  are  suited 
for  small  targets,  these  hold  promise  for  detecting  and  locating  breast  tumors  in  early  stages  of 
growth,  which  is  crucially  important  for  effective  treatment.  Further  work  involving  ex  vivo 
(model)  and  in  vivo  imaging  of  cancerous  breast  will  be  needed  to  establish  the  full  potential  of 
these  approaches. 
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Abstract:  A  time  reversal  optical  tomography  (TROT)  method  for  near- 
infrared  (NIR)  diffuse  optical  imaging  of  targets  embedded  in  a  highly 
scattering  turbid  medium  is  presented.  TROT  combines  the  basic  symmetry 
of  time  reversal  invariance  and  subspace-based  signal  processing  for 
retrieval  of  target  location.  The  efficacy  of  TROT  is  tested  using  simulated 
data  and  data  obtained  from  NIR  imaging  experiments  on  absorptive  and 
scattering  targets  embedded  in  Intralipid-20%  suspension  in  water,  as  turbid 
medium.  The  results  demonstrate  the  potential  of  TROT  for  detecting  and 
locating  small  targets  in  a  turbid  medium,  such  as,  breast  tumors  in  early 
stages  of  growth. 
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1.  Introduction 

Optical  imaging  of  targets  embedded  in  a  highly  scattering  turbid  medium,  such  as,  a  tumor  in 
a  breast,  is  a  challenging  problem  because  light  is  strongly  absorbed  and  scattered  by  the 
medium  leading  to  poor  signal-to-noise  ratio,  as  well  as,  loss  of  phase  coherence  and 
polarization.  As  a  consequence  distinct,  sharp  image  of  the  targets  may  not  be  formed  directly. 
Various  frequency-domain,  time -resolved,  and  steady-state  inverse  image  reconstruction  (HR) 
[1-5]  approaches  are  being  pursued  to  form  tomographic  images  using  diffusively  scattered 
light  measured  at  the  sample  boundary.  HR  is  an  ill-posed  problem  and  the  development  of 
reliable  and  fast  approaches  remains  a  formidable  task.  Recent  HR  algorithms,  such  as 
Newton-Raphson-Marquardt  algorithms  [6]  and  direct  linear  inversion  of  3-D  matrices  [7], 
are  time  consuming.  The  iterative  methods  [7,8]  may  not  ensure  that  the  obtained  result 
arrives  at  a  “global  minimum”  or  converges  to  a  “local  minimum”.  Still  the  potential  for 
developing  non-invasive  imaging  approaches  with  diagnostic  ability  motivates  the  ongoing 
diffuse  optical  tomography  (DOT)  research  using  NIR  light. 

Many  applications  require  rather  accurate  determination  of  location  of  target(s)  in  three 
dimensions.  For  example,  a  recent  study  involving  35,319  patients  underscores  the  influence 
of  primary  tumor  location  on  breast  cancer  prognosis  [9],  and  makes  it  imperative  that  DOT 
for  breast  cancer  detection  be  able  to  obtain  three-dimensional  (3-D)  location  of  the  tumor. 
While  two-dimensional  (2 -D)  HR  approaches  may  provide  only  lateral  positions,  3-D  HR 
approaches  attempt  to  retrieve  all  three  position  coordinates  of  the  target(s).  Various 
frequency-domain,  time-domain,  and  steady-state  DOT  approaches  have  addressed  the  target 
localization  problem  with  different  measures  of  success  [1-8].  Several  groups  have  paid 
particular  attention  to  retrieving  target  location.  Kepshire  et  al  developed  a  subsurface  DOT 
approach  to  obtain  location  information  of  absorbing  and  fluorescent  targets,  but  observed  the 
sensitivity  to  vary  nonlinearly  with  depth  [10].  Mohajerani  et  al  reported  a  fluorescent 
tomography  method  for  locating  fluorescent  targets  embedded  in  a  heterogeneous  medium 
using  partitioning  of  the  fluorophore  distribution  into  an  object  subspace  and  a  background 
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subspace  [11].  Godavarty  et  al  developed  another  fluorescent  tomography  approach  that  used 
a  hemispherical  breast  phantom,  near-infrared  light-induced  fluorescence  from  a  contrast 
agent,  and  finite  element  method-based  reconstruction  algorithms  to  obtain  target  information 
up  to  a  depth  of  2  cm  from  breast  phantom  surface  [12].  Zhao  et  al  introduced  a  layer-based 
sigmoid  adjustment  method  to  improve  depth  resolution  of  DOT  and  achieved  positioning 
error  within  3  mm  for  depths  from  10  to  30  mm  [13].  Optical  tomography  using  independent 
component  analysis  (OPTICA)  approach  developed  by  Xu  et  al  uses  multi-source  probing 
and  multi-detector  signal  acquisition  scheme  and  a  numerical  algorithm  based  on  independent 
component  analysis  of  information  theory  to  obtain  3-D  position  information  of  absorbing, 
scattering  and  fluorescent  targets  embedded  in  highly  scattering  turbid  media,  and  “model 
breast”  assembled  using  ex  vivo  human  breast  tissue  [14-17].  Co-registration  approaches  that 
use  another  modality,  such  as,  ultrasound,  magnetic  resonance  imaging,  and  x-ray 
mammography  for  locating  suspect  areas  and  DOT  for  obtaining  images  have  also  been 
introduced  [18-21]. 

In  this  article  we  report  on  the  development  of  a  time  reversal  optical  tomography  (TROT) 
[22-24]  approach  for  NIR  optical  imaging  of  target(s)  in  a  turbid  medium,  and  present  initial 
results  of  its  efficacy  using  both  simulated  and  experimental  data. 

Time  reversal  (TR)  invariance,  the  basic  symmetry  that  commonly  holds  in  microscopic 
physics,  forms  the  basis  for  macroscopic  TR  imaging.  TR  imaging  using  the  so-called  “time- 
reversal  mirrors”  (TRMs)  has  been  used  as  an  experimental  tool  in  acoustics  with  practical 
applications  in  medicine,  underwater  imaging,  and  nondestructive  testing  [25-28].  The 
theoretical  and  numerical  techniques  involved  in  time  reversal  have  been  used  for  applications 
involving  both  acoustic  waves  and  electromagnetic  waves  (radar)  [28-33]. 

Devaney  and  associates  developed  a  theoretical  framework  for  a  TR  imaging  method  with 
Multiple  Signal  Classification  (MUSIC)  for  finding  the  location  of  scattering  targets  whose 
size  is  smaller  than  the  wavelength  of  acoustic  waves  or  electromagnetic  waves  (radar)  used 
for  probing  the  homogeneous  or  inhomogeneous  background  medium  in  which  the  targets 
were  embedded  [34,35].  While  their  initial  focus  was  on  back-propagation  geometry  that  used 
coincident  acoustic  or  electromagnetic  transceiver  array  for  interrogating  the  targets,  they  later 
extended  the  formalism  to  transmission  geometry  where  sources  and  detectors  were  distinct 
and  separated  [36].  They  also  generalized  the  theory  which  was  based  on  distorted  wave  Born 
approximation  (DWBA)  to  account  for  multiple  scattering  between  the  targets  [37].  In  its 
basic  form  TR-MUSIC  found  target  location  from  knowledge  of  the  response  matrix  K,  which 
was  constructed  from  multi-static  data  collected  by  the  transceiver  array  [34,35].  TR-MUSIC 
provided  higher  spatial  resolution  than  the  conventional  TR  imaging,  especially  in  the  case 
where  targets  were  not  well  resolved  [34,35,38]. 

We  are  adapting  and  extending  the  TR-MUSIC  approach  to  the  optical  domain,  i.e.  to 
diffusive  optical  imaging  for  detecting  and  locating  targets  embedded  in  a  turbid  medium.  In 
this  paper,  TROT  is  studied  in  details  using  both  simulated  data  and  data  from 
transillumination  NIR  imaging  experiments  in  slab  geometry.  A  TR  matrix  is  obtained  by 
multiplying  the  response  matrix  formed  using  experimental  or  simulated  data  to  its  conjugate 
matrix.  The  leading  non-zero  eigenvalues  of  the  Hermitian  TR  matrix  determine  the  signal 
subspace  due  to  presence  of  the  targets.  The  signal  subspace  is  separated  from  the  noise 
subspace  using  an  Z-curve  method  [5,39,40].  The  vector  subspace  method,  MUSIC,  along 
with  Green’s  functions  calculated  from  an  appropriate  forward  model  for  light  propagation 
through  the  turbid  medium  is  then  used  to  determine  the  locations  of  the  targets.  The  MUSIC 
algorithm  judges  if  the  calculated  Green’s  function  vector  corresponding  to  a  location  in  the 
sample  is  mapped  into  the  signal  subspace  or  the  noise  subspace. 

Several  salient  features  make  TROT  attractive  and  potentially  more  promising  than  other 
HR  methods.  First  the  size  of  the  TR  matrix  is  much  smaller  than  those  used  in  other  HR 
approaches,  which  makes  solution  of  the  eigenvalue  problem  easier  and  faster.  Second,  to 
determine  locations  of  targets,  TR-MUSIC  approach  runs  the  program  over  all  voxels  only 
once,  and  there  is  no  need  to  carry  out  an  iterative  procedure  done  by  other  inverse 
approaches.  Other  HR  approaches  seek  to  determine  the  absorption  and  scattering  parameters 
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at  all  voxels  into  which  the  sample  is  divided.  The  process  is  iterative,  computationally 
intensive,  and  leads  to  a  solution  of  the  inverse  problem  that  is  not  unique  because  the 
problem  is  ill-posed,  even  when  there  is  no  noise.  In  contrast,  TROT  seeks  to  determine  the 
locations  of  the  targets  first  and  thereafter  retrieve  other  information,  such  as,  the  size  and 
optical  properties  of  the  limited  number  of  targets  in  the  medium,  which  requires  significantly 
less  computation  time.  The  focus  of  this  paper  is  on  finding  the  locations  of  targets. 

Our  result  using  simulated  data  shows  that  without  the  presence  of  noise  TROT 
determines  the  locations  of  the  embedded  targets  accurately  with  high  resolution.  TROT 
exhibits  promise  to  locate  targets  both  in  simulations  and  experiments  even  when  substantial 
noise  is  present.  Images  of  small  targets  obtained  by  this  approach  are  sharper  than  that 
obtained  by  other  HR  approaches. 

This  paper  is  organized  as  follows.  In  section  2,  the  formalism  of  the  TROT  approach  is 
presented.  In  section  3,  the  numerical  algorithm  of  TROT  is  described.  In  section  4,  the 
efficacy  of  the  formalism  is  tested  using  simulated  data.  Section  5  presents  the  results  when 
the  formalism  is  applied  to  experimental  data  using  Intralipid-20%  suspension  in  water  as  the 
highly  scattering  turbid  medium.  Section  6  discusses  the  results. 

2.  Formalism 

2. 1  Diffusion  approximation,  perturbation  method  and  response  matrix 

The  starting  point  for  the  TROT  formalism  is  the  diffusion  approximation  [41-43]  of  the 
radiative  transfer  equation  (RTE)  [44,45].  The  perturbation  in  the  light  intensity  distribution 
due  to  small  inhomogeneities  (targets)  embedded  in  a  homogeneous  medium,  to  the  first  order 
Bom  approximation,  can  be  written  as  [46,47] 

A<t>(rd,rs )  =  ~iG(rd,r)Sua(r)cG(r,rs)d3r 

-\SD(r)cVrG(rd,r)-VrG(r,rs)d2r,  (1) 

where  rs,  rd,  and  r  are  the  positions  of  a  point-like  source  of  unit  power,  detector  and  target, 
respectively;  G(r,rs )  and  G(rd,r )  are  the  Green’s  functions  that  describe  light  propagations 

from  the  source  to  the  target  and  from  the  target  to  the  detector,  respectively;  Spa  is  the 
difference  in  absorption  coefficient  and  SD  is  the  difference  in  diffusion  coefficient  between 
the  targets  and  the  background  medium;  and  c  is  the  light  speed  in  the  medium. 

A  multi-source  interrogation  and  multi-detector  signal  acquisition  scheme  is  used  to 
acquire  transillumination  data,  from  which  the  difference  in  the  light  intensity  distribution  due 
to  the  targets,  =  is  found,  where  $  is  the  light  intensity  distribution  measured  on 

the  sample  boundary  with  targets  embedded  in  the  scattering  medium  and  is  ideally  the 
light  intensity  distribution  without  the  targets,  which  in  practice  is  approximated  by  an 
“average”  over  all  the  multi-source  measurements.  A  response  matrix  K  is  constructed  with 
-A(j)  ,  to  describe  the  transport  of  light  from  different  sources  through  the  embedded  objects 
to  the  array  of  detectors  [22,36]. 

For  small,  point-like  absorptive  targets,  the  matrix  elements  can  be  rewritten  in  a  discrete 
form  as: 


1V1 

K,,  =  X°d  (nXm)rmGs (Xm,rj),i  =  1,2,- -,Nd;j  =  1,2, 


,N. 


(2) 


where  rm  =  Sjua(Xm)cSVm  is  the  optical  absorption  strength  of  the  mth  target,  SVm  is  the 
volume  of  mth  target,  rh  rd  and  Xm  are  locations  of  the  ith  detector,  jth  source  and  mth  target, 
respectively.  Due  to  the  reciprocity  of  light  propagation  in  the  medium,  G(r,r')  =  G(r\r)  . 
Thus, 
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(3) 


and 


1V1 

KtJ  =YGd  (Xm,rt)Tm<?  (Xm,rj), 


M 

K  =  {K^  =  J^d(Xm)rmgsT  (Xm), 

m=\ 


(4) 


where  g5(r)  and  gd(r)  are  Green’s  function  vectors  (GFVs)  associated  with  the  source  array 
and  detector  array,  respectively.  GFVs  are  defined  as 


gs,  (r)  =  [(?  (r, ,r),Gs  (r2,r),-,G*  (rv  ,r)f ,  (5a) 

gd  (r)  =  \Gd  (, r„r),Gd  (r2,r),-- -,Gd  (rNj,r)f,  (5b) 


where  the  superscript  T  denotes  transpose;  and  Ns,  Nd  and  M  are  the  numbers  of  sources, 
detectors  and  targets,  respectively.  It  is  assumed  the  number  of  targets  is  less  than  the  number 
of  sources  and  detectors,  M  <  min( Nd ,  Ns )  .  It  also  holds  that  KT  =  {K  .. }  describes  light 

propagation  from  the  positions  of  detectors  through  the  medium  and  targets  to  sources. 

For  a  homogeneous  background  medium,  the  rank  R  of  matrix  K,  is  equal  to  the  dimension 
of  the  source  array  vector  space  Qs  spanned  by  gs  (rm ) ,  and  also  equal  to  the  dimension  of  the 

detector  array  vector  space  Qd  spanned  by  gd  (rm )  ,  where  Qs  ^  CNs  and  Qd  cC^  .  For 
absorptive  targets,  R  is  equal  to  the  number  of  targets  M. 

Similar  forms  of  the  response  matrix  and  GFVs  can  be  obtained  for  scattering  targets.  As 
the  dot  product  in  the  second  term  of  Eq.  (1)  implies,  each  scattering  target  is  represented  by 
three  components  coexisting  at  one  location.  The  elements  of  the  K  matrix  for  L  scattering 
target  may  be  written  as 


L 


1=1 


=  Yj,  Z  da<7'(rt,Xl)daG’(Xl,rJ),  (6) 

l-\  a={x,  y,z} 


where  rl  =SD(Xl)cSVl  is  the  optical  scattering  strength  of  the  Ith  target.  The  K  matrix  for 
scattering  targets  can  be  written  in  a  manner  similar  to  that  for  absorptive  targets: 

K  =  t  Z  ^(X^^g/iX,).  (7) 

1=1  a={x,y,z } 


The  Green’s  function  for  a  slab  geometry  is  [16,47] 


-I  oo 

G(r,r')-G(r'r)-—  I 


(  -Krt  -la \  ^ 

e  e 


V  rk  rk  j 


K  =  [(-x_^')2  +{y~y')2 +(z+z'+  2M)2] 


(8a) 

(8b) 


where  k  =  [_(jua  -  ia>  /  c)  /  D^1  in  frequency  domain  with  amplitude  modulation  frequency  co, 
and  k  =  0,±  1,  ±2,  •  •  - .  The  extrapolated  boundaries  of  the  slab  are  located  at  z  =  0  and 
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z  =  d  =  L  +  2ze  ,  respectively,  where  L  is  the  physical  thickness  of  the  slab  and  the 

extrapolation  length  ze  is  determined  from  the  boundary  condition  of  the  slab  [48,49]. 

Under  ideal  conditions,  when  all  three  scattering  components  of  each  of  the  L  scattering 
targets  are  well-resolved,  the  rank  of  K  contributed  by  L  scattering  targets  is  3 L.  In  practice, 
four  components  (one  for  absorption  and  three  for  scattering)  are  calculated  for  each  target, 
since  the  targets  may  have  both  scattering  and  absorptive  characteristics,  or  the  exact  nature 
may  not  be  known  a  priori.  The  dominant  characteristic  is  used  to  label  the  target  as 
absorptive  or  scattering  in  nature. 

2.2  Point  Spread  Functions 

If  light  emitted  by  a  source  of  unit  power  at  target  position  X  propagates  in  the  sample 
medium,  the  signal  measured  by  the  detector  array  at  the  sample  boundary  is  Gd(rf,X)  .  The 
signal  is  then  “time -reversed”  and  back-propagated  with  the  Green’s  function  of  the 
background  medium.  TR  operation  is  phase  conjugation  in  Fourier  domain  [28,50].  So  the 
signal  evaluated  at  r  is  [34] 

Hd  (r,X)  =  f ]Gd  (r,r )  Gd*  (r,X)  =  g/  (r ) g/  (X ) 

i= 1 

=  gPx)gd{r)  =  {gd{X),gd{r)),  (9) 

where  *  denotes  phase  conjugate,  j*  denotes  adjoint,  and  <  •  >  denotes  inner  product. 
Hd(r,X)  is  the  detector  array  point  spread  function  (PSF).  A  source  array  PSF  can  be 
similarly  formed  as 

Hs  (r,X)  =  g,  (X  f  gs  (r)  =  (gs  (X),gs  (r)).  (10) 

Due  to  the  time  reversal  assumption,  Hd(r,X)  peaks  at  r  —  X  ,  so  it  can  be  considered  as 
an  image  of  the  source  at  X  formed  by  the  TR  detector  array.  PSF  vanishes  when  r  is  far  away 
from  X.  A  similar  interpretation  can  be  used  for  Hs(r ,  X)  . 

2.3  Time  reversal  and  MUSIC 

The  TR  matrix  may  be  constructed  to  represent  light  propagation  from  sources  to  detectors 
and  back  denoted  by  TSDDS,  or  to  represent  light  propagation  from  detector  positions  to  source 
positions  and  back  denoted  by  TDSSD,  a  consequence  of  the  reciprocity  of  light  propagation 
[29,34,38,50,51].  For  frequency-domain  data,  TSDDS=KtK  ,  and  TDSSD  =  (KT )? KT  =  K  KT  , 
where  response  data  matrix  K  is  formed  using  modulated  intensities,  instead  of  the  field  with 
phase  information  used  in  the  conventional  TR.  For  CW  measurements,  TSDDS  =  KT K  ,  and 
TDSsd  -  KKT  (K  is  real  and  only  includes  intensity  values). 

Since  TSDDS  and  TDSSD  are  Hermitian  (  =  T ),  they  have  complete  sets  of  orthonormal 

eigenvectors  v7  (  j  =  1,---,NS)  and  ut  (i  =  1, — ,  Nd  ),  with  a  common  set  of  non-negative  real 
eigenvalues.  For  M  <  min( Ns ,  Nd  )  absorptive  targets  without  the  presence  of  noise,  the  rank 
of  TSdds  and  Tdssd  is  M.  The  eigenvalues  Aj  >  0  ,  when  /  =  1,---,M  ,  and  /l  .  «  0  ,  when 
j  =  M  +  l,---,Ns  for  Tsdds  and  /  =  M  + 1, •  •  • , Nd  for  TDSSD-  The  eigen  system  {y. , ud ,  Aj  >  0}  , 
/  =  1,-  •  -  ,M  ,  is  related  to  the  targets.  The  TR  matrix  TSDDS  can  be  written  as  [22,34] 

M  M 

tsdds =TZC*AUxMMm))S;Um)gPx,).  (H) 

m=lm'=l 

Subsequent  formalism  may  be  different  depending  on  whether  the  targets  are  “well 
resolved”  or  “poorly  resolved.” 
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2.3.1  Well-resolved  targets 

If  the  mth  and  m'th  targets  (  )  are  well  resolved,  defined  by  the  conditions: 

Hs(Xm,Xm,) « 0  and  Hd(Xm,Xm,) « 0  ,  z.e.  the  GFVs  at  Xm  and  are  orthogonal, 

<  gd{Xm),gd(Xm)  >=  Hd(Xm„XJ  =  \\gd  (Xm)f  Smm, .  So  we  have 

M 

4*»=ZkJ2||&(^)||  (i2) 


where  ||*||  denotes  L2  norm  [52].  The  eigenvectors  of  TSDDS  are  proportional  to  the  phase 
conjugate  of  the  GFVs  associated  with  the  M  targets  [22,34],  i.e. 


TsnnsS :  (Xm )  =  |rj2  \\gd  (xrn  )f  \\gs  (Xm  )f  g/  (XJ. 
The  eigenvectors  are 


v,-  = 


Ss 

M 

Ss 

M 

(13) 


(14) 


with  eigenvalues  |^(^y)|  |&(^7-)|  »  /  =  1,---,M.  Thus  TSDDS  is  a  projection 

operator  that  projects  a  vector  onto  the  conjugate  of  the  source  array  vector  space  Qs .  The  jth 
non-zero  eigenvalue  Xj  is  directly  related  to  the  optical  strength  r7  of  the  jth  target.  Similar 
equations  can  be  derived  for  TDSSD,  which  is  a  projection  operator  for  the  conjugate  of  the 
detector  array  vector  space  Qd  .  The  eigenvectors  of  T ossd  are 


Sd 


(15) 


/  =  1,*  •  -  ,M  ,  with  the  same  eigenvalues  as  TSdds- 

Therefore,  for  well-resolved  targets,  the  target  locations  can  be  determined  by  the  inner 
product  [22,34,36,51] 


gx(X 


*.(*,) 


(16a) 


or 


Vi 


’•sAx,) 


J 


& 


Ml' 


(16b) 


where  Xp  is  a  test  target  position,  which  is  the  position  of  any  voxel  in  the  sample  space,  y/f 
and  if/jd  peak  when  Xp  is  the  position  of  the  jth  target.  In  the  classical  TR  imaging 
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[25,29,38,50],  for  ideally  resolved  targets,  each  eigenvector  of  the  TR  operator  can  be  used  to 
focus  on  one  particular  target.  Here  y/f  and  y/f  represent  focusing  of  signals  from  the  source 
and  detector  planes  on  to  the  target  position,  respectively.  Use  of  the  eigenvectors  Vj  and  itj, 
/  =  1,-  •  -  ,M  ,  ensures  that  the  jth  target  is  sorted  out.  When  TSDDS  and  source  array  vector  space 
( Tdssd  and  detector  array  vector  space)  are  used,  we  call  the  scheme  SDDS  (DSSD).  Both 
source  and  detector  arrays  can  be  considered  simultaneously  to  locate  the  target  by  calculating 


V,  =v,dv?  = 


UuOINMu) 


(17) 


/  =  1,*  •  -  ,M  ,  which  is  computationally  equivalent  to  a  process  that  light  emitted  from  a  virtual 
source  of  unit  power  at  a  test  target  position  Xp,  propagates  to  the  TR  source  array  and  back  to 
a  true  target  position  Xj\  then  it  is  re-emitted  and  further  propagates  to  the  TR  detector  array 
and  back  to  the  original  position  Xp.  y/7  peaks  when  the  test  target  position  Xp  coincides  with 
the  true  target  position  Xj  associated  with  the  jth  eigenvector. 

2.3.2  Poorly-resolved  Targets  and  MUSIC 

When  the  targets  are  too  close  to  each  other  or  the  sources  and/or  detectors  are  significantly 
sparse,  the  targets  are  considered  to  be  poorly  resolved  and  the  GFVs  at  Xm  and  Xm<  are  not 
orthogonal.  In  such  cases,  the  eigenvectors  v7  and  Uj  do  not  correspond  one-to-one  with  the 
GFVs  associated  with  target  positions  Xm  (  /,/w  =  1,---,M).  The  image  resolution  degrades 
because  of  contributions  from  multiple  targets.  To  solve  this  problem,  the  subspace-based 
method,  MUSIC  was  implemented  with  TR  [34,36,51].  MUSIC  algorithm  is  based  on  the  idea 
that  although  the  vectors  characterizing  the  targets  are  no  longer  orthogonal  with  each  other, 
they  are  all  located  in  the  signal  subspace,  which  is  orthogonal  to  the  noise  subspace. 

The  orthonormal  sets  {v*}  (  f  =  !,-••, Ns)  and  {u*}  (  j  =  l,--,Nd)  span  the  spaces  CNs 
and  CNd  associated  with  the  source  and  detector  arrays,  respectively.  While  {v  *}  and  {u*}  , 
with  Aj  >  0  ,  form  the  signal  subspaces  on  the  source  and  detector  arrays,  Ss  ={v*}  and 
Sd  =  {u* }  (  /  =  1,---,M  ),  respectively;  {v*}  and  {a*}  ,  with  A.  « 0  ,  form  the  noise 
subspaces,  A fs  =  {v* }  (  j  =  M  +  1,---,NS)  and  Nd  =  {u* }  (  f  =  M  +  l,---,Nd),  respectively. 

Thus  CNs  =  Ss  ®  J\fs  and  CNd  -  Sd  ©  J\fd  [36,51].  Since  the  dimensions  of  the  signal 
subspaces  Ss  and  Sd  and  of  the  GFV  spaces  Qs  and  Qd  are  all  equal  to  M,  Qs  =  Ss  and 
Qd  =  Sd  [51].  The  GFVs,  gs(Xm)  and  gd(Xm),  m  =  1  ,  are  linear  combinations  of  Vj* 

and  uj*,  /  =  1,---,M,  respectively.  Therefore,  gs(Xm)<=Ss  and  gd{Xm)^Sd  ,  w  =  l,---,M, 
associated  with  mth  target  are  orthogonal  to  v  *  g  A fs  (  /  =  M  + 1,  •  •  • ,  Ns  )  and  u*  g  A fd 
(  /  =  M  + 1,  •  •  • ,  Nd ),  respectively: 

(v/,&  (^))  =  v/&  (xm)  *  0,  j=M  + 1,-  •  -,N„  (18a) 

(«/>&  (^J)  = =«/&  (V„)  »0,  j  =  M  +  \,.  ■  -,Nd.  (18b) 

The  locations  of  targets  can  be  determined  by  calculating  the  following  squared  sum  of 
inner  products: 

Qs(*P)=  £  <19a) 

j=M +1 
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(19b) 


&(*,)=  I  Wsd{xPi- 

j=M+ 1 

g5(X^)  and  vanish  when  the  test  target  position  Xp  is  a  true  target  position.  Similar 

to  Eq.  (17),  Q  =  QsQd  can  be  calculated  with  both  source  and  detector  arrays  considered 
simultaneously.  An  alternative  approach  to  accentuate  a  target  position  is  to  plot  a  pseudo 
spectrum  defined  as 

Ps(Xp)  =  \\gs(Xp)f  /\Qs(Xp)\  (20a) 

associated  with  the  source  array,  or 

pd{xP)=h{xPi  I\qMp)\  (20b) 

associated  with  the  detector  array,  or 

p(xr)-pXxr)pAxr)  <20c> 

associated  with  both  source  and  detector  arrays  [22,34,36,51],  where  and 

lk„(*,)||2  are  used  for  normalization.  The  poles  of  the  pseudo  spectrum  correspond  to  target 

locations.  These  MUSIC  pseudo  spectra  can  also  be  used  to  locate  well -resolved  targets. 

Since  the  dimension  of  the  signal  subspace  is  generally  much  smaller  than  that  of  the  noise 
subspace,  it  is  preferred  that  in  Eq.  (19)  and  Eq.  (20),  the  signal  subspace  is  used  rather  than 
the  noise  subspace  for  ease  of  computation.  Using  the  properties  of  the  projection  operators 
associated  with  the  source  and  detector  arrays  [22,34,36,51],  Qs(Xp)  and  Qd{X  )  can  be 
calculated  as 

Qs{Xp)  =  \\gs{Xp)f-±\v/gs(Xp)\2,  (21a) 

7=1 

ii  m2  ^  i  i2 

a(u)=Mu)|  I  •  (21b) 

7=1 

When  the  targets  are  embedded  in  a  non-uniform  medium,  or  when  there  is  significant 
noise  present,  the  noise  or  false  targets  contribute  significantly  to  the  eigenvalues.  The  near¬ 
zero  and  non-zero  eigenvalues  are  not  as  well  separated  as  when  there  are  no  noise.  In  this 
case,  the  rank  of  the  TR  matrix  is  larger  than  the  number  of  targets  M.  The  TR  matrix  may 
even  be  full  rank.  However,  as  long  as  M  is  less  than  min (Ns,Nd)  and  eigenvalues 
contributed  by  the  noise  and  false  targets  are  smaller  than  those  contributed  by  the  real  targets 
with  a  threshold  e  ,  the  target  positions  can  be  obtained  using  a  pseudo  spectrum  [36,51] 
associated  with  the  source  array, 

(22) 

II  ||2  I  T  |2 

where  Qs{Xp)^.<e  =  ~  A^\vj  Ss(^P)\  •  Pseudo  spectra  associated  with  the  detector 

Aj  >£ 

array  or  with  both  source  and  detector  arrays  can  also  be  obtained  similarly.  In  practice,  the 
threshold  is  selected  to  separate  the  signal  and  noise  subspaces  using  a  method  similar  to  L- 
curve  regularization  [39]. 


#153647  -  $15.00  USD  Received  31  Aug  2011;  revised  29  Sep  2011;  accepted  1  Oct  2011;  published  21  Oct  2011 
(C)  201 1  OSA  24  October  201 1  /  Vol.  19,  No.  22  /  OPTICS  EXPRESS  21965 


When  scattering  targets  are  concerned,  the  GFVs  dag  (a  =  x,y,z),  associated  with  the 
test  target  position  Xp  will  be  used  to  calculate  the  pseudo  spectrum.  For  a  target  with  both 
absorption  and  scattering  properties  at  the  wavelength  of  probing  light,  one  GFV 
corresponding  to  absorption  constructed  as  g  and  three  GFVs  corresponding  to  scattering 
target  constructed  with  dag ,  (a  =  x,y,z),  are  used  to  calculate  the  pseudo  spectrum  over 
every  voxel.  Ideally,  for  an  absorptive  and  scattering  target  four  pseudo-values  will  be 
obtained  for  every  target  position.  If  the  dominant  value  corresponds  to  the  absorptive  (any  of 
the  scattering)  GFV  the  target  will  be  identified  as  absorptive  (scattering)  in  nature. 

3.  Algorithm 

Implementation  of  TROT  to  locate  targets  embedded  in  a  highly  scattering  turbid  medium 
involves  the  steps  outlined  below.  For  simplicity,  the  sizes  of  source  array  and  detector  array 
are  assumed  to  be  the  same,  i.e.,  Nd  =  Ns  =  N  . 

(a)  A  response  matrix  K  with  size  N  x  N  is  constructed  using  experimental  data  (or 

estimated  data  in  simulation).  Data  consist  of  the  perturbations  in  the  light  intensity 
distribution  due  to  the  targets,  A(/>  =  $  -  ,  where  (f)  is  the  light  intensity  distribution 

measured  on  the  sample  boundary  with  targets  embedded  in  the  scattering  medium 
and  c/)0  is  ideally  the  light  intensity  distribution  without  the  targets.  In  practice,  (/)0  is 
approximated  by  an  “average”  over  all  the  multi-source  measurements,  while  in 
simulation  it  can  be  estimated  without  such  approximation. 

(b)  A  detector  array  TR  matrix,  TDSSD  =  KKT  with  size  N  x  N  for  CW  measurements  is 
constructed.  All  the  eigenvalues  and  the  eigenvectors  of  the  TDSSD  matrix  are 
computed.  The  eigenvectors  are  orthogonal  to  each  other.  It  is  to  be  noted  that  in  this 
procedure  we  only  deal  with  a  matrix  of  dimension  N ,  not  a  matrix  of  dimension  of  N 
x  N  as  done  in  traditional  inverse  procedures. 

(c)  The  non-zero  eigenvalues  of  TDSSD  belonging  to  the  signal  subspace  are  separated 
from  the  near-zero  eigenvalues  belonging  to  the  noise  subspace  using  the  T-curve 
method  [5,39,40]. 

(d)  MUSIC  approach  [34,36,51]  is  next  used  to  determine  the  locations  of  the  targets  as 

follows,  (i)  The  3-D  medium  is  divided  into  a  certain  number  of  voxels.  A  detector 
array  GFV,  gd  ( Xp  ) ,  associated  with  an  absorptive  test  target  position  Xp  at  pth  voxel 
is  calculated  using  Diffusion  Approximation  of  RTE.  Other  proper  forward  models 
could  be  used  as  well.  In  order  to  check  if  gd(Xp )  is  located  in  the  signal  subspace 

or  in  the  noise  subspace,  a  pseudo  spectrum  associated  with  the  detector  array  is 
computed  using  Eq.  (20b),  where  M  is  the  dimension  of  the  signal  subspace  found  in 
step  (c).  If  gd(X  )  is  located  in  the  signal  subspace,  the  corresponding  pseudo  value 

P(Xp )  in  Eq.  (20b)  will  become  a  maximum,  (ii)  Pseudo  spectra  are  also  calculated 

using  the  other  three  GFVs,  dagd  (X  ) ,  (  a  =  x, y,z )  for  scattering  property,  (iii)  All 

pseudo  values  are  put  together  and  sorted  in  a  descending  order.  Since  the  leading 
pseudo  values  at  Xp  are  associated  with  targets  and  specific  GFVs,  the  positions  of 
the  embedded  targets  and  their  nature  (absorptive  or  scattering)  are  determined.  The 
pseudo  spectrum  in  the  whole  sample  space  can  be  used  to  plot  pseudo  tomographic 
images. 

In  this  approach,  only  a  single  run  is  needed  for  calculating  the  pseudo  spectrum  and  no 
iterative  procedure  is  involved,  which  makes  it  faster  and  computationally  less  intensive  than 
the  traditional  HR  approaches.  Similar  procedure  can  be  used  for  application  of  TROT  when 
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Nd  *  Ns .  The  pseudo  spectrum  associated  with  either  the  source  array,  or  the  detector  array, 
or  both  source  and  detector  arrays,  as  outlined  in  Eq.  (20)  can  be  used  to  obtain  target 
positions. 

It  is  instructive  to  compare  the  computational  complexity  of  TROT  formalism  with  that  of 
typical  iterative  methods.  For  a  typical  iterative  method,  an  equation  b  =  Wx  is  solved  to  find 
the  inhomogeneities  (targets),  where  W  is  a  weight  matrix  with  size  NdNs  x  N,  N  is  the  number 
of  voxels,  b  is  an  NdNs  x  1  vector  describing  the  perturbation  in  the  detected  light  intensity  due 
to  the  presence  of  inhomogeneities,  and  x  is  the  perturbation  or  variation  in  the  optical 
properties  from  the  background  values  with  dimension  of  N  x  1.  The  computational 
complexity  is  typically  0(NdNsN2)  for  a  single  iteration.  The  computational  complexity  of 
TROT  is  much  smaller  than  that  for  even  one  iteration  of  an  iterative  method.  For  the  SDDS 
scheme,  the  complexity  for  TROT  is  0(NdNs2)  if  NdNs  >  NNk,  and  0(NsNNk)  otherwise,  where 
Nk  is  the  dimension  of  the  signal  subspace.  In  the  DSSD  scheme,  the  complexity  is  0{NsNd)  if 
NdNs  >  NNk,  and  0(NdNNk)  otherwise.  TROT  does  not  involve  any  iteration. 

In  the  following  sections,  TROT  will  be  tested  using  simulation  and  experimental  data. 

4.  Testing  TROT  with  simulated  data 

To  test  the  efficacy  of  the  TROT  approach,  we  first  consider  a  rather  challenging  task  of 
detecting  and  locating  six  targets  embedded  in  a  simulated  sample  which  is  a  40 -mm  thick 
uniform  scattering  slab.  Its  absorption  and  diffusion  coefficients  are  /ia  =  1/300  mm-1  and  D  = 
1/3  mm,  respectively.  The  incident  CW  beam  was  step-scanned  in  an  x-y  array  of  41  x  41  grid 
points  with  a  step  size  of  2  mm  on  the  input  plane  covering  an  80  mm  x  80  mm  area.  Fight 
transmitted  from  the  opposite  side  (output  plane)  was  recorded  at  41  x  41  grid  points  covering 
the  same  area.  No  random  noise  was  added. 

The  six  (M=  6)  point-like  absorptive  targets,  with  absorption  coefficient  difference  of  Ajua 
=  0.01  mm-1  from  the  background,  were  placed  at  A  (24  mm,  26  mm,  9  mm),  B  (38  mm,  38 
mm,  15  mm),  C  (38  mm,  38  mm,  21  mm),  D  (40  mm,  38  mm,  21  mm),  E  (44  mm,  42  mm,  21 
mm)  and  F  (30  mm,  30  mm,  31  mm),  respectively.  The  origin  (0  mm,  0  mm,  0  mm)  was 
located  at  the  upper-left  corner  of  the  input  boundary  (source  plane)  of  the  sample.  The 
medium  was  divided  into  40  x  40  x  20  voxels,  with  each  voxel  of  size  2  mm  x  2  mm  x  2  mm. 
As  can  be  seen  from  the  assigned  coordinates,  targets  C  and  D  are  located  at  two  adjacent 
voxels,  and  are  close  to  target  E,  and  these  three  targets  are  located  in  the  same  z  layer. 
Consequently,  targets  C  and  D  are  expected  to  be  very  difficult  to  resolve,  and  hard  to 
distinguish  from  target  E.  Target  B  and  C  have  the  same  lateral  position  x  and  y,  and  different 
depths.  Target  A  is  close  to  the  source  plane,  while  F  is  close  to  the  detector  plane. 

Using  the  Diffusion  Approximation  of  the  RTE  as  the  model  for  light  propagation  in  slab 
geometry,  signals  arising  from  light  propagation  from  the  source  array  to  the  detector  array 
through  medium  with  and  without  the  targets  were  calculated.  The  difference  between  the  two 
sets,  which  is  the  perturbation  due  to  the  targets,  was  used  as  the  “simulated  data”.  The  size  of 
the  K  matrix  is  N  x  N  =  1681  x  1681.  The  TR  matrix  T=  KKr  was  constructed.  Then,  1681 
eigenvalues  and  1681  eigenvectors  of  T  were  found. 

The  first  seven  (7)  computed  eigenvalues  in  a  descending  order  of  magnitude  are  listed  in 
the  first  column  of  Table  1.  The  leading  twenty  eigenvalues  are  plotted  in  Fig.  1(a)  on  a 
logarithmic  scale.  The  first  six  (6)  eigenvalues  are  at  least  10  orders -of-magnitude  higher  than 
the  7-th  and  other  smaller  eigenvalues.  Hence,  the  dimension  of  the  signal  subspace  and  the 
number  of  targets  are  determined  to  be  six.  The  pseudo  spectrum  (consisting  of  40  x  40  x  20 
x  4  elements)  was  calculated  using  the  M  eigenvectors  in  the  signal  subspace.  The  values  of 
elements  in  the  pseudo  spectrum  were  sorted  in  a  descending  order.  The  seven  leading  pseudo 
values  are  listed  in  Table  1  with  the  corresponding  positions  of  voxels.  The  six  peaks  are 
found  to  be  associated  with  the  GFVs  for  absorptive  targets.  Namely,  the  corresponding  six 
targets  are  identified  to  be  absorptive  targets. 

All  six  large  pseudo-values  are  located  at  the  exact  known  target  locations  and  their  values 
are  approximately  9  orders -of-magnitude  larger  than  those  at  their  neighborhood  locations.  A 
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2-D  slice  of  the  pseudo  spectrum  on  z  =  21  mm  plane  is  shown  in  Fig.  1(b),  showing  the 


locations  of  the  three  difficult  targets. 


(a) 


Fig.  1.  (a)  A  plot  of  first  twenty  (20)  eigenvalues  on  logarithmic  scale,  (b)  2-D  slice  of  the 
pseudo  spectrum  on  z  =  21  mm  plane  showing  the  location  of  the  three  difficult  targets 
described  in  the  text.  Similar  2-D  slices  were  also  obtained  for  z  =  9-mm,  15-mm,  and  31 -mm 
planes  (not  shown). 


Table  1.  Eigenvalues,  pseudo  spectrum  and  the  corresponding  positions 


Leading 

Eigenvalues 

Poles  of  Pseudo 
Spectrum 

Retrieved  Position 

(x,  y,  z)  mm 

Known  Position 

(x,  y,  z )  mm 

2.6697E-010 

1.591  IE  +  015 

(44,  42,21) 

(44,  42,21) 

1.1722E-01 1 

8.6376E  +  014 

(38,38,  15) 

(38,38,  15) 

4.0081E-013 

7.9559E  +  014 

(38,38,21) 

(38,38,21) 

3.6676E-014 

7.2328E  +  014 

(40,38,21) 

(40,38,21) 

5.2629E-016 

6.3010E  +  014 

(24,  26,  9) 

(24,  26,  9) 

6.4837E-017 

2.1159E  +  014 

(30,  30,31) 

(30,  30,31) 

2.8337E-039 

2.4353E  +  005 

(38,  38,  19) 

With  the  highly  encouraging  result  from  simulation  even  for  a  considerably  challenging 
task,  we  proceeded  to  test  the  approach  for  the  realistic  situation  of  detecting  and  locating 
targets  from  experimental  data. 

5.  Testing  TROT  using  Experimental  Data 

5. 1  Experimental  materials  and  methods 

Three  different  experiments  with  three  different  samples  were  carried  out  to  test  the  efficacy 
of  the  TROT  approach  in  detecting  and  locating  targets  in  a  turbid  medium.  All  three  samples 
used  a  250  mm  x  250  mm  x  60  mm  transparent  plastic  container  filled  with  Intralipid-20% 
suspension  in  water  as  the  background  medium.  The  concentration  of  Intralipid-20%  was 
adjusted  to  provide  an  estimated  [53,54]  absorption  coefficient  pa  -0.003  mm-1  at  790  nm, 
and  a  transport  mean  free  path  lt  ~  1  mm,  which  were  similar  to  the  average  values  of  those 
parameters  for  human  breast  tissue,  while  the  cell  thickness  of  60  mm  was  comparable  to 
thickness  of  a  typical  compressed  breast. 

In  the  first  experiment,  the  depth  (position  along  z-axis)  of  an  absorptive  target  was  varied 
to  explore  how  the  accuracy  of  position  estimate  depended  on  depth.  The  target  was  a  glass 
sphere  of  diameter  ~9  mm  filled  with  ink  dissolved  in  Intralipid-20%  suspension  in  water  (jus 
was  adjusted  to  be  the  same  as  that  of  the  background  medium,  while  jua  «0.013  mm-1  was 
about  3  times  higher  than  that  of  background  medium). 
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In  the  second  experiment,  the  separation  between  two  absorptive  targets  was  varied  to  test 
how  close  those  could  be  and  yet  be  resolved  as  separate  objects.  Both  the  targets  were  similar 
to  the  target  in  the  first  experiment. 

In  the  third  experiment,  the  depth  of  a  scattering  target  was  varied  to  explore  the  efficacy 
of  TROT  in  locating  and  characterizing  a  scattering  target.  The  target  was  a  glass  sphere  of 
diameter  10  mm  filled  with  Intralipid-20%  suspension  in  water  to  provide  a  transport  mean 
free  path  lt  of  0.25  mm,  and  a  scattering  coefficient  jus&  11  mm-1. 

A  multi-source  interrogation  and  multi-detector  signal  acquisition  scheme,  shown  in  Fig. 
2,  was  used  to  acquire  data.  A  100-mW  790-nm  diode  laser  beam  was  used  to  illuminate  the 
samples.  A  1024  x  1024  pixels  charge  coupled  device  (CCD)  camera  equipped  with  a  60-mm 
focal-length  camera  lens  was  used  on  the  opposite  side  of  the  sample  to  detect  the  transmitted 
light  on  the  boundaries  of  the  slab  samples  (detector  plane).  The  pixel  size  was  24  pm.  The 
multi-source  illumination  scheme  was  realized  by  scanning  the  sample  across  the  laser  beam 
in  a  two-dimensional  x-y  array  of  grid  points  using  a  computer-controlled  translation  stage. 
The  first  and  third  samples  were  scanned  across  the  laser  beam  in  an  array  of  9  x  9  grid 
points,  and  the  second  sample  was  scanned  in  an  array  of  1 5  x  1 1  grid  points,  with  a  step  size 
of  5  mm  in  all  cases.  The  scanning  and  data  acquisition  processes  were  controlled  by  a 
personal  computer  (PC).  Raw  transillumination  images  of  the  sample  were  recorded  by  the  PC 
for  each  scan  position,  and  stored  for  subsequent  analysis.  A  typical  image,  which  is  a  2-D 
intensity  distribution,  is  shown  in  the  right  side  of  Fig.  2. 


Fig.  2.  A  schematic  diagram  of  the  experimental  arrangement  for  imaging  objects  embedded  in 
a  turbid  medium.  (Key:  CCD  =  charge  coupled  device,  PC  =  personal  computer)  Inset  (below) 
shows  the  2-D  array  in  the  input  plane  that  was  scanned  across  the  incident  laser  beam,  and 
inset  (right)  shows  a  typical  raw  image. 

While  we  have  scanned  the  sample  and  kept  the  source  fixed  in  the  experiments  reported 
here,  a  more  clinically  relevant  approach  would  be  to  scan  the  source  and  fix  the  sample.  In 
the  experimental  arrangement,  the  source  scanning  may  be  accomplished  by:  (a)  delivering 
the  beam  using  an  optical  fiber,  and  translating  the  delivery  end  of  the  fiber  in  an  x-y  array 
using  a  computer-controlled  translation  stage;  or  (b)  raster  scanning  the  laser  beam  using  two 
orthogonal  (x-y)  galvanometers.  The  main  change  in  the  processing  of  data  would  involve 
alignment  of  the  images  so  that  laser  beam  positions  are  overlapped  before  averaging  to 
generate  the  background  image. 

5.2  Data  Processing  and  Analysis 

From  each  image,  a  region  of  interest  was  cropped  out  and  then  every  5x5  pixels  in  the 
cropped  image  were  binned  to  one  pixel  to  enhance  the  signal-to-noise  ratio.  The  background 
image  was  generated  by  taking  an  average  of  the  original  images  for  all  scan  positions,  which 
is  a  reasonable  approximation  since  for  most  of  the  scan  positions  the  target(s)  is  (are)  not 
along  the  direction  of  the  incident  beam.  Then  the  background  image  was  also  cropped  and 
binned  corresponding  to  the  region  of  interest  for  each  scan  position.  Perturbation  in  the  light 
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intensity  distribution  A(j)  due  to  targets  in  each  image  was  found  by  subtracting  the 
background  image  from  each  individual  image.  The  response  matrix  was  constructed  using 
the  light  intensity  perturbations,  -A</>  .  The  TR  matrix  was  generated  by  multiplying  the 

response  matrix  by  its  transpose  for  our  continuous -wave  (CW)  probing  scheme.  The 
eigenvalue  equation  was  solved  and  the  signal  subspace  was  selected  and  separated  from  the 
noise  subspace.  MUSIC  was  then  used  to  calculate  the  pseudo  spectrum  for  all  voxels  in  the 
3-D  space  of  the  sample.  For  each  voxel,  four  pseudo  values,  one  for  absorption  and  three  for 
scattering  as  described  in  Algorithm  step  (d),  were  calculated.  The  voxel  size  was  0.77  mm  x 
0.77  mm  x  1  mm.  By  sorting  the  pseudo  spectrum  in  a  descending  order,  the  target(s)  were 
located. 

The  voxel  size  to  be  used  in  reconstruction  and  its  relation  to  the  target  size  is  an  important 
consideration.  In  general,  smaller  voxels  provide  reconstruction  of  higher  resolution  at  the 
cost  of  increased  computational  time.  Finer  details  of  an  extended  target  may  be  obtained 
using  smaller  voxels.  Decreasing  the  voxel  size  indefinitely  may  not  improve  resolution 
because  of  the  diffusive  nature  of  light  propagation  in  the  turbid  medium.  However,  the 
computation  time  increases  dramatically,  which  has  been  observed  by  other  researchers  [55]. 
The  optimal  voxel  size  for  a  given  reconstruction  problem  will  depend  on  factors,  such  as, 
target  size,  experimental  geometry,  and  noise  level. 

Since  the  signal  used  in  image  reconstruction  is  taken  to  be  the  difference  between  the 
image  recorded  for  a  scan  position  and  the  background  image,  estimation  of  the  background 
image  is  an  important  issue.  This  is  a  common  problem  for  every  diffuse  optical  imaging 
modality  using  perturbation  method,  and  needs  further  elaboration.  We  accumulated  data  in 
the  transillumination  slab  geometry,  and  generated  the  background  image  by  averaging 
images  for  all  scan  positions  after  proper  alignment  with  respect  to  the  incident  source 
position.  This  averaging  method  for  generating  background  image  worked  well  for  small 
targets  that  we  used  in  our  experiments,  as  the  ratio  of  sample  volume  to  target  volume  was 
quite  high  (-500:1).  This  volume  ratio  for  breast  and  a  tumor  in  early  stages  of  growth  will 
also  be  substantially  high  for  the  averaging  method  to  be  applicable.  Assuming  a  scenario 
where  the  volume  ratio  is  substantially  smaller  than  in  above  examples,  a  modified  approach 
would  be  to  select  recorded  images  which  were  minimally  affected  by  embedded  targets  for 
averaging  [56].  As  long  as  the  targets  only  occupy  a  limited  volume  within  the  host  medium,  a 
clean  background  image  can  be  generated  in  this  fashion.  It  should  also  be  noted  that  while 
estimation  of  target  optical  properties,  such  as  absorption  coefficient  and  scattering 
coefficient,  are  sensitively  dependent  on  background  image  estimation,  estimation  of  target 
positions  are  not  so  sensitive. 

Several  alternative  ways  of  generating  background  image  are  suggested  in  the  literature. 
One  experimental  approach  is  to  record  image  using  a  phantom  that  has  the  same  average 
optical  properties  as  the  sample,  such  as  human  breast  [57].  Along  the  same  line,  image  of  the 
healthy  contralateral  breast  taken  under  the  same  experimental  conditions  as  that  of  the 
suspect  breast  may  be  used  as  background  image  for  breast  imaging  [58].  Some  authors  have 
suggested  acquiring  data  at  a  wavelength  for  which  the  target(s)  and  the  background  have 
identical  optical  properties  for  assessing  the  background,  e.g.,  measurement  using  805-nm 
light  for  which  hemoglobin  and  oxyhemoglobin  have  the  same  absorption  coefficient  may 
serve  as  background  to  image  hemoglobin  oxygenation  [59].  Still  another  approach  is  to 
compute  the  background  using  an  appropriate  forward  model  [18].  Any  of  these  approaches 
may  be  employed  for  generating  the  background  image  for  use  with  the  TROT  formalism 
presented  here. 

The  geometries  commonly  used  in  DOT  include  slab,  cylindrical,  hemispherical,  and 
semi-infinite;  and  different  source -detector  combinations  have  been  used  to  record  images  in 
these  geometries.  As  long  as  multiple  source -detector  combinations  provide  multiple  angular 
views  of  the  sample  the  TROT  formalism  can  be  adapted  to  obtain  target  location  for  these 
geometries.  TR  imaging  and  TR-MUSIC  was  originally  developed  for  reflection 
(backscattering)  geometry  that  used  coincident  transceiver  array  to  detect  the  return  signal 
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[28-30,34,35].  With  requisite  modification  in  the  experimental  arrangement  TROT  would  be 
suited  for  use  in  the  reflection  geometry. 

5.3  Results 

5.3.1  Single  absorptive  target  at  different  depths 

In  the  first  experiment,  only  one  target  was  used,  the  lateral  (x,  y)  position  of  the  target  was 
kept  the  same  at  (25.5  mm,  24.7  mm),  while  seven  different  depths  (position  along  z-axis)  of 
15  mm,  20  mm,  25  mm,  30  mm,  35  mm,  40  mm  and  45  mm  were  used.  The  eigenvalue 
spectrum  plotted  on  logarithmic  scale  for  the  target  at  z  =  30  mm  is  shown  in  Fig.  3.  Similar 
eigenvalue  spectra  were  obtained  for  other  cases. 


Fig.  3.  A  semi-log  plot  of  eigenvalue  spectrum  with  first  40  leading  eigenvalues  for  the  target 
at  z  =  30  mm. 

Both  SDDS  and  DSSD  pseudo  spectra  were  calculated  using  Eq.  (20).  The  target  was 
identified  as  an  absorptive  target.  In  the  DSSD  pseudo  spectrum,  the  absorptive  pseudo  value 
at  the  peak  position  is  ~41  times  of  the  scattering  pseudo  value  associated  with  dzgd  ,  and 
even  larger  than  those  associated  with  dxgd  ,  and  dygd  ,  as  shown  in  Table  2.  Similarly  in  the 
SDDS  scheme,  the  absorptive  pseudo  value  at  the  peak  position  is  -33  times  of  the  scattering 
pseudo  value  with  dzgs ,  and  much  larger  than  the  other  two. 

Table  2.  Pseudo  values  associated  with  absorptive  and  scattering  components  at  the  peak 

position 


Scheme  with 

Absorptive 

Scattering  pseudo  value 

GFVfe) 

pseudo  value 

8yS  dzg 

DSSD  (gd) 

1305.0 

1.0 

1.0  31.7 

SDDS  fe) 

2729.3 

14.0 

1.1  81.6 

Three-dimensional  tomographic  images  were  generated  using  the  whole  absorption  pseudo 
spectrum  for  all  voxel  positions  in  the  sample.  The  left  pane  of  Fig.  4(a)  shows  a  tomographic 
image  for  the  target  at  z  =  30  mm.  The  spatial  profiles  in  the  x,  y  and  z  directions,  shown  in  the 
right  pane  of  Fig.  4(a)  were  used  to  assess  the  target  location.  Similar  images  were  generated 
for  other  depths.  The  retrieved  target  positions  are  compared  with  known  positions  in  Table  3. 

As  is  evident  from  Table  3,  when  DSSD  scheme  was  used,  the  TROT-assessed  lateral 
positions  (x,  y)  were  within  0.6  mm  of  the  known  values,  which  is  an  excellent  agreement. 
The  accuracy  of  the  z-position  was  found  to  be  optimal  when  the  target  was  located  in  the 
middle  plane  of  the  sample,  and  deteriorated  when  the  target  was  closer  to  the  source  plane  or 
the  detection  plane.  When  using  SDDS  scheme,  the  TROT-assessed  lateral  positions  were 
also  within  0.6  mm  of  the  known  positions,  except  for  z  =  40  mm  and  45  mm,  where  the  error 
in  y  direction  was  1 .2  mm  and  2  mm,  respectively.  However,  remarkable  improvement  in  the 
accuracy  of  the  z-position  estimation  was  observed,  the  error  Az  being  within  0.5  mm  for  all 
cases  except  for  z  =  35  mm,  where  the  error  was  1.5  mm.  We  ascribe  the  superior 
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performance  of  the  scheme  using  TSDDs->  to  the  much  larger  size  of  the  detector  array  (1024  x 
1 024)  than  that  of  the  source  array  (9  x  9)  used  in  the  experimental  arrangement. 

Table  3.  Positions  of  one  target  located  at  different  depths 


Known  Positions 
v,  y,  z  (mm) 

DSSD  Scheme 

SDDS  Scheme 

Retrieved 
Positions 
v,  y,  z  (mm) 

Error 

Ax,  Ay,  Az 
(mm) 

Retrieved 
Positions 
v,  y,  z  (mm) 

Error 

Ax,  Ay,  Az 
(mm) 

25.5,24.7,  15 

24.9,  24.4,  17.5 

0.6,  0.3,  2.5 

24.9,  25.2,  15.5 

0.6,  0.5,  0.5 

25.5,  24.7,  20 

25.7,  24.4,21.5 

0.2,  0.3,  1.5 

24.9,  25.2,  20.5 

0.6,  0.5,  0.5 

25.5,  24.7,  25 

25.7,  24.4,  26.5 

0.2,  0.3,  1.5 

25.7,  24.4,  24.5 

0.2,  0.3,  0.5 

25.5,  24.7,  30 

25.7,  24.4,  30.5 

0.2,  0.3,  0.5 

25.7,  25.2,  29.5 

0.2,  0.5,  0.5 

25.5,24.7,35 

25.7,  25.2,33.5 

0.2,  0.5,  1.5 

24.9,  24.4,  36.5 

0.6,  0.3,  1.5 

25.5,  24.7,  40 

24.9,  25.2,  36.5 

0.6,  0.5,  3.5 

24.9,  25.9,  40.5 

0.6,  1.2,  0.5 

25.5,  24.7,  45 

24.9,  25.2,  39.5 

0.6,  0.5,  5.5 

24.9,  26.7,  45.5 

0.6,  2.0,  0.5 

(b)  10  20  30  40  50 

x(mm) 


10  20  30  40  50 
x(mm) 


1200 

300 

400 


x  10 

14 

■3 

-2 


14 


1200 

BOO 

400 


Fig.  4.  Pseudo  image  of  the  target  (left  pane)  and  corresponding  spatial  intensity  profiles  (right 
pane)  when  the  target  is  located  at  z  =  30  mm:  (a)  experimental  data;  (b)  simulation  without 
any  added  noise;  and  (c)  simulation  with  20%  Gaussian  noise  added.  The  pseudo  values  are 
calculated  using  Eq.  (20). 

It  should  be  noted  that  the  choice  of  either  DSSD  or  SDDS  scheme  depends  on 
experimental  parameters,  such  as,  the  number  and  density  of  sources  and  detectors,  and  does 
not  depend  on  the  characteristics  of  the  background  medium.  When  more  detectors  than 
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sources  are  used  and  inter-detector  spacing  is  small,  SDDS  would  provide  better  resolution 
than  DSSD,  and  vice  versa.  However,  due  to  the  diffusive  nature  of  light  propagation  in  the 
turbid  medium,  increasing  the  numbers  and  decreasing  the  spacing  of  the  sources/detectors 
beyond  a  limit  may  not  always  improve  the  results. 

While  the  target  position  could  be  obtained  from  the  experimental  data,  it  was  observed 
that  the  difference  between  the  smaller  eigenvalues  in  the  signal  subspace  and  the  larger 
eigenvalues  in  the  noise  subspace  were  not  as  pronounced  as  observed  in  the  simulation  in 
Section  4.  To  assess  the  effect  of  noise  and  to  what  extent  noise  may  be  present  in  the 
experimental  data;  we  generated  simulated  data  mimicking  the  experimental  conditions,  and 
added  different  noise  levels.  The  lateral  positions  were  (25.5  mm,  24.7  mm)  and  all  seven  z- 
positions  (depth)  of  15  mm,  20  mm,  25  mm,  30  mm,  35  mm,  40  mm,  and  45  mm  were  tested. 
Typical  pseudo  images  generated  for  z  =  30  mm  without  and  with  20%  Gaussian 
multiplicative  noise  to  compare  with  the  experimental  result  are  shown  in  Fig.  4(b)  and  Fig. 
4(c),  respectively.  Simulated  data  provided  the  known  position  coordinates. 

The  simulated  spatial  profiles  with  zero  added  noise  are  much  sharper  than  the  profiles 
obtained  from  experimental  data,  or  from  simulated  data  with  20%  added  Gaussian  noise. 
Broadening  of  spatial  profile  is  an  indication  of  the  uncertainty  in  determination  of  position 
coordinates.  Results  from  simulation  show  that  uncertainty  in  position  determination  increases 
with  added  noise,  and  that  experimental  data  behave  in  a  way  similar  to  simulated  data  with 
added  noise. 

5.3.2  Resolving  two  absorptive  targets 

In  the  second  experiment  using  two  targets  the  depth  (z)  and  height  (y)  were  kept  same  (z  =  30 
mm,  y  =  26.0  mm),  while  three  different  center-to -center  separations,  Ax  of  -12.6  mm,  17.6 
mm,  and  27.6  mm,  between  them  along  the  x-axis  were  considered.  A  cross-sectional  pseudo 
image  of  the  targets  when  separated  by  a  center-to-center  distance  of  27.6  mm,  generated 
using  the  pseudo  spectrum  is  shown  in  the  left  pane  of  Fig.  5(a).  Figure  5(b)  shows  a  similar 
image  for  the  separation  12.6  mm  (separation  between  nearest  edges  -4  mm).  A  similar  image 
for  the  separation  17.6  mm  was  also  obtained  (not  shown  in  the  figure).  The  profiles  in  the  x,  y 
and  z  directions  through  the  right  target  are  shown  in  the  right  panes  of  Fig.  5(a)-Fig.  5(c). 
These  profiles  were  used  to  assess  locations  of  the  targets,  and  the  separation  between  the  two 
targets.  In  all  cases,  the  targets  were  determined  to  be  absorptive,  because  peaks  occurred  in 
the  pseudo  spectrum  with  the  GFVs  corresponding  to  absorption  property. 

The  known  and  retrieved  positions  from  the  experiments  and  separations  Ax  between  the 
two  targets  appear  in  Table  4.  In  all  the  cases,  the  two  targets  were  resolved,  even  when  their 
center-to-center  separation  was  12.6  mm  apart,  nearest  sides  separated  by  only  -4  mm.  For  all 
retrieved  positions,  the  maximum  error  in  the  lateral  positions  is  3.0  mm,  and  the  maximum 
error  in  the  axial  positions  is  1.5  mm.  The  errors  in  the  lateral  positions  increase  as  the  targets 
get  closer.  We  ascribe  this  increase  in  error  in  the  lateral  position  to  the  crosstalk  between  the 
two  targets,  the  peak  due  to  one  target  being  affected  by  the  other.  The  shift  in  the  peaks  is 
also  affected  by  noise.  When  the  two  targets  are  very  close  or  significant  noise  is  present,  the 
two  peaks  merge,  so  that  the  two  targets  are  not  resolved.  This  behavior  was  confirmed  in 
simulations. 

The  results  were  compared  with  simulated  data  using  similar  conditions.  For  the  more 
challenging  case  of  two  targets  located  at  z  =  30  mm  and  separated  by  12.6  mm,  exact  target 
locations  were  found  when  no  noise  was  added.  With  10%  noise  present,  the  positions  of  the 
two  targets  were  found  to  be  (39.0  mm,  24.8  mm,  29.0  mm)  and  (30.0  mm,  24.8  mm,  29.0 
mm)  with  target  separation  9.0  mm,  compared  to  12.6  mm  (known)  and  6.9  mm  retrieved 
from  experiment.  The  pseudo  image  and  the  corresponding  profiles  through  the  right  target 
are  shown  in  Fig.  5(c).  Similar  images  were  also  obtained  for  the  left  target.  The  retrieved 
separation  between  the  two  targets  in  simulation  with  1 0%  noise  was  smaller  than  the  actual 
separation.  But  the  error  was  less  than  the  experimental  result.  However,  when  20%  noise  was 
added,  the  two  peaks  merged  (not  shown  here). 
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Fig.  5.  (a)  (Experiment):  TROT  generated  cross-section  pseudo  image  when  the  targets  are 
separated  by  27.6  mm  is  shown  in  the  left  pane  and  pseudo-value  profiles  through  the  right 
target  along  x,  y  and  z  directions  are  shown  in  the  right  pane,  (b)  (Experiment):  TROT 
generated  cross-section  pseudo  image  when  the  targets  are  separated  by  12.6  mm  is  shown  in 
the  left  pane  and  the  corresponding  spatial  profiles  through  the  right  target  along  x,  y  and  z 
directions  are  shown  in  the  right  pane,  (c)  (Simulation):  TROT  generated  cross-section  pseudo 
image  when  two  targets  are  separated  by  12.6  mm  is  shown  in  the  left  pane  and  the 
corresponding  pseudo-value  profiles  are  plot  in  the  right  pane.  In  simulation  10%  Gaussian 
noise  is  added  for  comparison  with  the  experimental  results.  P  is  pseudo  value  calculated  using 
Eq.  (20). 


Table  4.  Positions  of  two  targets  separated  with  different  distances 


Known 

Target  # 

Known 

Retrieved 

Error 

Retrieved 

Separation 
[Ax  (mm)] 

Position 
[x,  y,  z  (mm)] 

Position 
[x,  y,  z  (mm)] 

(mm) 

Separation 
[Ax  (mm)] 

12.6 

1 

2 

27.6,  26.0,  30 
40.2,  26.0,  30 

30.3,24.4,28.5 
37.2,  25.2,  29.5 

2.7,  1.6,  1.5 
3.0,  0.8,  0.5 

6.9 

17.6 

1 

25.1,26.0,  30 

26.4,  24.4,  28.5 

1.3,  1.6,  1.5 

14.6 

2 

42.7,  26.0,  30 

41.0,  25.2,  29.5 

1.7,  0.8,  0.5 

27.6 

1 

20.1,26.0,  30 

19.5,  25.2,  29.5 

0.6,  0.8,  0.5 

27.6 

2 

47.7,  26.0,  30 

47.1,25.2,  30.5 

0.6,  0.8,  0.5 

The  limits  on  the  size  of  targets,  separation  between  the  targets,  and  the  target-to- 
background  contrast  ratio  that  are  needed  to  detect  and  locate  the  targets  depend  on  noise 
level,  experimental  parameters  (such  as,  number  and  concentration  of  sources  and  detectors), 
and  ultimately  on  the  diffuse  nature  of  light  propagation  in  the  turbid  medium.  Coordinated 
experimental  work  and  numerical  modeling  will  be  needed  to  assess  those  limits. 
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5.3.3  Single  scattering  target  at  different  depths 


x.y.z  (mm) 


Fig.  6.  Pseudo  image  of  the  target  (left  pane)  and  corresponding  spatial  intensity  profiles  (right 
pane)  when  the  target  is  located  at  z  =  30  mm:  (a)  experimental  data;  (b)  simulation  with  20% 

Gaussian  noise  added.  P  is  pseudo  value  calculated  using  Eq.  (20). 

The  experiment  involving  the  third  sample  is  the  same  as  the  first  one  except  that  the  target 
was  scattering  in  nature.  The  scattering  target  was  a  10-mm  diameter  glass  sphere  filled  with 
Intralipid-20%  suspension  in  water,  whose  concentration  was  adjusted  to  provide  lt  =  -0.25 
mm,  fis=  11.3  mm-1.  The  same  scanning  and  data  acquisition  scheme  was  used  as  for  the 
absorptive  targets  and  the  following  z-positions  of  the  target  were  used:  15  mm,  20  mm,  25 
mm,  30  mm,  35  mm,  40  mm,  and  45  mm.  DSSD  scheme  was  used  to  calculate  the  pseudo 
spectrum.  A  cross-section  pseudo  image  and  the  corresponding  spatial  profiles  are  displayed 
in  Fig.  6(a)  when  z  =  30  mm.  It  is  compared  to  the  simulation  results  with  20%  Gaussian  noise 
(Fig.  6(b)).  The  lateral  (x,  y)  spatial  profiles  of  the  pseudo  image  generated  using  simulated 
data  are  considerably  wider,  while  the  axial  (z)  spatial  profile  is  narrower  than  those  obtained 
using  experimental  data,  and  the  peak  values  from  the  two  cases  are  of  the  same  order.  The 
retrieved  target  positions  are  listed  in  Table  5.  SDDS  scheme  was  also  used  and  provided  with 
similar  results. 


Table  5.  Positions  of  one  scattering  target  located  at  different  depths 


Retrieved  Positions  Error 

[x,  y,  z  (mm)]  [Ax,  Ay,  Az  (mm)] 

0.8,  1.4,  3.5 


Known  Positions 
[x,  y,  z  (mm)] 


25.7,  24.5,  15 

25.7,  24.5,  20 

25.7,  24.5,  25 

25.7,  24.5,  30 

25.7,  24.5,  35 

25.7,  24.5,  40 

25.7,  24.5,  45 


24.9,  25.9,  18.5 
27.2,  26.7,  20.5 

25.7,  26.7,  23.5 

24.9,  25.2,  32.5 

24.9,  25.2,  36.5 

24.9,  25.9,41.5 

24.9,  25.9,  45.5 


1.5,  2.2,  0.5 
0.0,  2.2,  1.5 
0.8,  0.7,  2.5 
0.8,  0.7,  1.5 
0.8,  1.4,  1.5 
0.8,  1.4,  0.5 


In  Fig.  5,  and  more  prominently  in  Fig.  6,  the  image  resolution  seems  better  for 
experimental  data  than  simulated  data.  Since  the  peak  values  and  bandwidth  of  lines  (the 
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poles)  in  the  pseudo  spectrum  depend  strongly  on  the  noise,  this  difference  in  image 
resolution  is  presumably  due  to  lower  noise  level  in  the  experiments  than  that  used  in 
simulations. 

A  comparison  of  experimental  results  for  scattering  and  absorptive  targets  validate  the 
common  notion  that  it  is  more  challenging  to  locate  and  image  scattering  targets  than 
absorptive  targets  in  a  highly  scattering  medium.  Also  the  lateral  (x,  y)  positions  are 
determined  with  higher  accuracy  than  the  axial  (z)  position.  Overall  the  TROT-retrieved  target 
positions  are  in  good  agreement  with  the  known  positions. 

6.  Discussion 

The  article  presents  the  development  of  time  reversal  imaging  approach  with  subspace 
classification,  MUSIC  in  the  optical  domain.  The  results  from  experiment  and  simulation 
show  that  TROT  is  a  faster  and  less  computation  intensive  approach  for  detecting  small 
targets  in  highly  scattering  turbid  media  and  determining  their  locations  in  3-D  than  other 
inverse  image  reconstruction  techniques.  While  the  dominant  features  in  the  pseudo  spectrum 
are  related  to  the  square  of  the  difference  between  the  absorption  (scattering)  coefficient  of  the 
targets  and  that  of  the  background,  the  approach  does  not  directly  determine  these  parameters. 
It  is  common  for  HR  approaches  to  estimate  the  optical  properties  of  every  voxel  in  the 
sample  and  identify  target(s)  from  differences  of  these  properties  between  the  sample  and  the 
target(s),  which  is  a  considerably  computation  intensive  undertaking.  On  the  contrary,  TROT 
identifies  the  targets  as  poles  of  the  pseudo  spectrum  and  focuses  on  determining  their 
positions,  which  do  not  require  as  much  computation  time.  Other  HR  approaches  involve 
iteration,  while  TROT  is  non-iterative.  In  TROT  the  data  dimension  is  lower  compared  to 
other  HR  approaches,  which  enables  analysis  and  utilization  of  very  large  data  sets.  These  two 
features  together  make  TROT  faster.  Fast  image  reconstruction  algorithms  are  of  particular 
interest. 

It  was  observed  that  lateral  (x,  y)  positions  are  better  determined  than  the  depth  (z).  Also 
the  spatial  profile  is  more  spread  out  along  z  compared  to  that  along  x,  y.  We  ascribe  this 
difference  to  fewer  data  along  z-direction  compared  to  those  along  x-y  planes.  Addition  of 
another  set  of  data  with  light  incident  and  signal  collected  perpendicular  to  the  z-direction  is 
expected  to  further  improve  resolution  in  this  dimension.  Even  without  that  addition,  TROT 
determines  the  target  position  well. 

While  we  have  used  slab  geometry  and  CW  illumination,  the  TROT  approach  may  be  used 
for  other  geometries  (such  as,  cylindrical,  and  spherical),  different  types  of  illumination  (e.g. 
frequency  domain  and  pulsed)  and  different  models  for  light  propagation  through  the  medium. 
Application  and  adaption  of  the  TROT  formalism  to  inhomogeneous  media  and  extended 
targets  may  require  careful  consideration  of  several  factors.  In  a  non-uniform,  inhomogeneous 
medium,  structures  other  than  the  desired  targets  may  appear  as  “false  targets”  and  may 
interfere  with  identification  of  “real  targets”.  However,  as  long  as  the  contributions  to  the 
signal  by  any  false  target  is  smaller  than  that  made  by  real  targets,  TROT  with  MUSIC  will  be 
useful  in  detecting  and  locating  targets,  by  choosing  a  proper  threshold  to  separate  the  signal 
and  noise  subspaces.  What  is  even  more  important,  expected  wavelength  dependence  of  the 
target  spectroscopic  properties  could  be  used  to  assess  the  difference  between  the  real  and 
false  targets  in  experiments  using  multi-wavelength  interrogation  of  the  sample. 

The  TROT  formalism  presented  in  this  article  is  particularly  suited  for  point-like  targets 
requiring  fewer  eigenvectors  in  the  signal  subspace  to  construct  a  pseudo  spectrum.  However, 
for  extended  finite-size  targets,  the  formalism  needs  to  be  modified  and  much  more 
eigenvectors  may  be  needed  to  calculate  the  pseudo  spectrum  [40,60,61].  These  interesting 
problems  for  further  study  are  currently  being  pursued. 
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Multi-functional  tumor-targeting  nanocomposites  and  time-reversal  optical  imaging  for  early 
detection  of  breast  cancer  and  prevention  of  micro-metastases 

Abstract  of  a  pre-proposal  submitted  to  the  Idea  Award  (collaborative  option)  ofBCRP  2011  by 
S.  K.  Gayen  (Initiating  PI)  and  Valeria  Balogh-Nair  (Collaborating  PI) 

City  College  of  New  York 

Research  Idea 

The  objective  of  the  proposed  research  is  to  develop  multi-functional,  tumor-targeting  nanocomposites 
and  near-infrared  (NIR)  optical  imaging  approaches  for  early  detection  of  breast  cancer,  and  prevention 
of  micro-metastases  that  are  respo  nsible  for  majority  of  breast  cancer  m  ortality.  The  nanocomposite 
synthesis  will  use  a  m  ultivalent  dendrimer  (a  class  of  organic  m  acromolecules)  platform  to  incorporate 
fluorescent  moieties  ( e.g .,  semiconductor  quantum  dots,  or  gold/silver  nanoparticles)  as  contrast  agents 
for  imaging,  and  chemokine  mimics  for  selective  targeting  of  cancer  cells  and  prevention  of  metastases. 
The  high  affinity  of  che  mokine  mimic’s  ligands  for  chemokine  CXCR4  receptors  will  ensure  selective 
delivery  of  the  nanoco  mposite  to  the  ta  rget.  A  multi-source  NIR  pr  obing  and  multi-de  tector  signal 
acquisition  arrangement  along  with  a  tim  e-reversal  image  reconstruction  algorithm  will  be  used  for  fast 
detection  of  tumor  and  determination  of  its  location  in  three  dimensions. 

Impact  and  Innovation 

The  major  impact  of  the  proposed  research  is  that  it  has  th  e  potential  to  provide  a  modality  for  early 
detection  of  breast  cancer,  and  for  prevention  of  metastases,  the  major  cause  of  mortality.  A  broader 
potential  impact  is  that  the  approach  could  be  extende  d  to  other  cancers  using  chem  okine  mimics 
directed  to  other  chemokine  receptors. 

The  project  is  highly  innovative  in  many  ways.  First,  the  synthesis  process  brings  together  the  concept  of 
multivalency,  the  idea  of  chemokine  mimics  for  prevention  of  cancer  cell  migration,  and  use  of 
dendrimers  as  stabilizers  and  nanoreactors  for  contrast  agent  synthesis.  A  combination  of  these  novel 
ideas  is  proposed,  to  the  best  of  our  knowledge  for  the  first  time,  as  a  multi-prong  approach  to  fight  the 
menace  of  breast  cancer.  Second,  the  affinity  of  the  chemokine  mimics  to  seven  helix  chemokine 
CXCR4  receptors  obviates  the  need  for  a  targeting  vector.  The  chemokine  mimics  will  play  the  dual 
role  of:  (1)  “homing  devices”  for  selective  delivery  of  the  nanocomposite  to  the  tumor  sites;  and  (2) 
“prevention  agents”  interacting  with  the  chemokine  receptor  sites  to  inhibit  metastases.  Third,  the 
efficacy  of  noninvasive  NIR  imaging  approach  for  early  detection  and  potential  diagnosis  will  be 
significantly  enhanced  through  design  of  efficient  contrast  agent.  Fourth,  the  idea  of  time-reversal 
optical  tomography  (TROT)  is  a  new  paradigm  in  diffusive  optical  tomographic  (DOT)  imaging.  While 
currently  pursued  DOT  approaches  are  iterative  and  computation  time  intensive,  TROT  is  non-iterative 
and  will  be  faster,  which  is  a  necessary  condition  for  real-time  imaging.  TROT  is  designed  for  detecting 
and  locating  small  targets  and  will  be  suited  for  early  detection  when  tumors  are  small.  Finally,  the 
dendrimer  platform  with  multivalent  surface  will  enable  incorporation  of  moieties  that  enhance  other 
imaging  modalities,  such  as,  magnetic  resonance  imaging,  enabling  development  of  sought-after  dual  or 
multimodal  imaging  modules. 


